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ABSTRACT MONTEREY CAL OVATE 


The problem of obtaining parametric models for linear 
and nonlinear systems based on observations of the input 
and output of the system is one of wide ranging interest. 
For linear systems, moving average (MA) and autoregressive 
(AR) models have received considerable attention and, based 
on the Levinson algorithm, a number of very powerful methods 
involving lattice filter structures have been developed to 
obtain the model solutions. For nonlinear systems the 
Volterra series model which is a nonlinear extension of the 
moving average model is frequently used. 

The purpose of this research is to extend these tech- 
niques to more general linear and nonlinear models. Using 
the equation error formulation, lattice solution methods 
in batch processing and adaptive form are developed for both 
Single and multichannel autoregressive moving average (ARMA) 
models for linear systems and Volterra series models for 
nonlinear systems. A nonlinear extension of the ARMA model 
is also considered and is shown in some cases to remedy 
problems encountered in Volterra modeling of nonlinear sys- 
tems. Lattice methods are also developed for the nonlinear 
ARMA model and it is shown that the methods obtained for 
linear ARMA modeling follow as a special case of the non- 
linear results. 

Experimental verification of the methods developed for 


Single channel linear ARMA modeling is presented and used 





to explore the characteristics of the lattice solution 
techniques. The results clearly indicate that the lattice 
Methods are extremely powerful, capable of producing highly 


accurate system models using short runs of data. 
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I. INTRODUCTION 


Traditionally, man has attempted to create models of 
portions of his environment for two principal reasons: 

1. To gain insight and understanding as to their 

per lonine, 

2, As a prelude to taking some action such as attempting 

to exercise control over them. 
The field of physics for instance, is replete with examples 
where men have created models to study and explain phenomenon 
from planetary motion to the motion and even origin of sub- 
atomic particles. In designing machines, engineers routinely 
rely on models of the components they use to describe how they 
will function, and to obtain the desired results in the final 
Product. Economics is another field where the use of models 
abounds for such purposes as identifying, forecasting or 
Eo uns to direct trends. 

The scope of the modeling problem is quite broad be- 
Ane with a decision on the type of model to be used, what 
physical quantities to measure, how to estimate the para- 
meters of the model from the measurement, and finally a veri- 
Beerıon of the model, In the chapters that follow, one 
facet of this problem, that of estimating model parameters, 
will be explored in detail for a number of linear and non- 


linear models, 





Ewe tHE PROBLEM STATEMENT 

Enespramary coneern er this work is the determination of 
discrete time models for both linear and nonlinear, time 
invariant systems from sampled observations of the system 
Mmpucs and Outputs, The general approach underlying all of 
the models examined here assumes that the system to be 


modeled is described by the equation 


y(k) = Figli CO Fg Ly Ck-212 JtF, buck) ,y(k-1)] ШЕТ? 


where Fig» Fg and Fao are functions of past and present 


values of their arguments, and u(k) and y(k) are the system 
Mut and output respectively. This is depicted in Figure 


(1.1). A possible method for modeling this type of system 


is to create a model of exactaly the same configuration with 


T3 


OE o 30° 


operate the system and model in parallel with the same input 


menctions and F E assume a form for these functions, 
and adjust the parameters of the model functions to minimize 
the mean square error (MSE) between the model output y(k) 

and the system output. The symbol "^" is used here to indi- 
cate that the signalis an estimate of y(k). This is depicted 


in Figure 1.2 and is often referred to as direct form modeling 


since the assumed topology of the system is directly copied 


Script notation will be used to refer to quantities 
associated with the system while nonscript notation will be 
used for their corresponding approximants in the model. 


ENS 












u(k) pae 
Figure 1.1. The assumed form for systems to be 
modeled. T represents a unit delay. 
SYSTEM Ди 
u(k) D e(k) 
y(k) 


MODEL 





Figure 1.2. A direct approach to system modeling. 
by the model. Here the model output 1s given by 


y(k) SF, plu(k)I+F, Ey (-1)1+F, 9 Lu 00 ,y Ck-1)] a» 


ДЕЛІ 





and the error signal is referred to as the "output error". 


As an example, if the system is linear an appropriate model 


1s 
N 
Dr l[luck)] = ) a(i)u(k-i) (1.34) 
n 
1-0 
N 
Fogty(k-1) J = ) | 5(i) y(Xk-i) ue B) 
= 
F, fuk) ,y(k-1)] = 0 (1.30) 
(for linear models, Fagl’ will be zero). In general 


however, the direct form approach will have serious diffi- 


lies if either Fag d ЕНЕ №] аге monzero since the 


30 
past values of yOO used in these functions are themselves 
dependent on the model parameters. A mimimum mean square 
output error approach results ina system of highly non- 
linear simultaneous equations which must be solved to 

obtain the model parameters. 

To avoid these difficulties, the equation error approach 
gut emn 2) to system modeling (which uses different 
model forms in the NA Sic and synthesis phases) will be 
Ended to the problem. The analysis model is depicted in 
Figure 1.3 and differs from the direct form model in that 
Poo and Er are functions of past and present values of the 
delayed system output rather than the analysis model output. 


Ша 





e(k) 


2 
ANALYSIS 
MODEL 





Figure 1.3. The equation error approach for 
system modeling. 


For each of the models studied, a general form for the three 
functions is assumed and the parameters of the model (coef- 
Enis Of the functions) are set to obtain a MMSE solution. 
ech case, the MSE cost function is a quadratic function 
of the model parameters (due to both the equation error 
formulation and to the form chosen for the functions) with 

a unique minimum and therefore the solution is given by a 
system of linear equations. The synthesis model is of the 
same form assumed for the system in Figure 1.1 and uses the 


nections Е E and Pag determined during the analysis 


10? 720 


phase. ) 


> 





As an alternative to the topology shown in Figure 1.3 
it will occasionally be convenient to consider the error 
signal e(k) as the output of the analysis model rather than 
mime, prediction v. This can be accomplished in one of two 


ways by defining 


FjLy Ck) ] IAE Pjgly(k-1)] Geta) 
or 


F4Lu GO ,y Ce] = y(k)-F, О у(к-1)] (ЛО) 


and reformulating the analysis model as shown in Figures 
l.4a or 1.4b. These model forms are often referred to as 
prediction error models since their outputs are the errors 
in predicting y(k) rather than the predictions themselves. 
There are, however, no substantive differences between the 
modeling approaches depicted in Figure 1.3, l.4a and 1.4b. 
The equation error formulation can be generalized to 
multiple input multiple output systems as well (henceforth 
referred to as multichannel systems) by considering u(k) and 
y(k) as vectors of Q: input signals and Q, output signals 


Е 


апа Pro? Fg» 30? Р, апа E, ВЕС Го functions. The 
prediction error signal e(k) becomes a Qgrvector of signals 
and the model parameters can be set to minimize the trace 
Шаде Prediction error covariance matrix P = с Те ое(ю 12. 
Such generalizations have been developed to a degree in the 


multichannel filtering literature and will be investigated 


mether here. 
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Figure 1,4. Prediction error model forms. 
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It is important to keep in mind however, that while the 
eEguation error formulation can be used to find a model solu- 
tion, it is an indirect method as opposed to the direct form 
method which minimizes the mean square value of output error. 
The direct form model has been modified to obtain the equa- 
tion error analysis model so that the parameters can be ob- 
tained via the solution of systems of linear equations. The 
price paid for this simplification in the model analysis 
problem is that additive noise on the measured system out- 


Duüt will introduce a bias in the model coefficient estimates. 


ЕЕ OVERVIEW 

Chapter II along with appendices A through F provide a 
unified review of the existing background theory on minimum 
mean square equation error modeling of linear systems. The 
moving average (MA) and autoregressive (AR) models are pre- 
sented and their relative merits compared. In Section II.C. 
the Levinson algorithm [Refs. 9, 10 and 27] for the AR and 
MA models is developed, greatly simplifying the solution 
process for these models. Section II.D, then shows that the 
Levinson algorithm defines the AR and MA models in terms of 
ШЫ се filter structures. 

These lattice structures have received widespread atten- 
tion and have led to a host of new developments in modeling, 
Spectral estimation, filter structures and adaptive filtering. 
Examination of the properties of these forms have suggested 


a number of new methods for calculating model coefficients 





that offer increased accuracy, and in some cases guarantee 
model stability even in the presence of correlation estimates 
obtained by averaging over short time intervals [Refs. 5, 
ENS and 361. Applied by Burg [Ref. 5] to spectral esti- 
mation, these methods allow the determination of power 
spectra via AR modeling from very short runs of data without 
any need for the use of a window function. [In finite pre- 
cision arithemetic implementations, the lattice structures 
have been shown by Markel and Gray [Ref. 33] to be less sensi- 
tive to roundoff noise and coefficient quantization than 
direct structures and have led to the design of other 
structures that offer improved performance over conventional 
parallel realizations. Griffiths has shown that these 


lattices can be implemented adaptively 





and that they offer the potential for more rapid convergence 
than conventional LMS adaptive filters. Recently Morf [Refs. 
36, 37 and 38] has also used these lattice structures to 
implement a recursively updated deterministic least squares 
adaptive scheme. It is readily apparent therefore, that the 
E inal work of Levinsan and the lattice structures that 
have evolved from it have had an important impact on the 
field of digital signal processing, 

Ша Section II.E.., the multichannel generalization of many 
of the single channel AR and MA modeling results is presen- 
ted. After a discussion of the basic multichannel AR and MA 
models [Refs. 26 and 45], the multichannel version of the 


Levinson algorithm originally developed by Whittle [Ref. 56], 


17 





and Wiggins and Robinson [Ref, 61] is presented. A new 
form for the models is introduced and used here however, 
to facilitate the application of these results later in 
Various other modeling problems. Multichannel lattice 
structures are then derived and from them alternative 
solution methods for the modeling problems are developed. 

Biaally, in section II.F. the LMS adaptive algorithm 
[Ref. 58] is reviewed and the adaptive implementations of 
the lattice structures due to Griffiths are presented. 

In Chapter III, the more general autoregressive moving 
average model is presented using the equation error formu- 
lation attributed to Kalman [Ref. 23]. After a brief 
discussion of the model, new model transition formulas 
are developed showing how the ARMA model is related to the 
Simpler and less general AR and MA models. System input 
Signal requirements for the ARMA modeling process are ex- 
plored and it is shown that the power spectrum of the input 
Signal can be considered as a frequency dependent weighting 
Function in the model optimization. Then the main result 
of the chapter is presented. With suitable assumptions, 

a recursive in order solution method for ARMA modeling 

(the (ntl)=st order solution is obtained from the n-th 

order solution) is obtained based on the Levinson algorithm for 
multichannel AR models. From this, lattice solution methods 
for the ARMA model are developed in both batch processing 

and adaptive form. (Batch processing here refers to assuming 


ergodicity and estimating correlations with time averaging). 


10 





A similar result has recently been presented by Morf [Refs. 
37 and 38] with the assumption of a white noise input sig- 
nal to the system. The results presented here follow from a 
different approach without the assumption of a white noise 
input. Experimental results are also presented verifying 
the methods and theory, and showing their advantages (and 
disadvantages) over conventional ARMA modeling methods. 

The programs used in these simulations are listed in Appen- 
ШЕСІ). In Section III.F., and Appendix б, it is shown that 
these single channel methods readily extend to the multi- 
channel ARMA model, and as one would expect, can be obtained 
as a special case. 

In Chapter IV two types of nonlinear models, the Volterra 
series model and the new nonlinear ARMA model recently pro- 
posed by Parker [Ref. 64], are considered. After a brief 
discussion of the Volterra model, it is shown that the 
solution can be obtained using multichannel MA lattice methods 
if the regular form of the Volterra kernels is used in place 
of the conventional symmetric form. Then the nonlinear ARMA 
Medel 1s presented in Section IV.B. and it is shown that for 
many systems, this model can remedy the problem of the large 
number of terms (ideally infinite) required by the Volterra 
model to represent the system in much the same way that the 
ARMA model solved the problem arising in the MA model. In 
section IV.B.2 it is also shown that by using the regular 
form, the solution for the nonlinear ARMA model can be 
obtained using multichannel AR lattice methods and that the 


linear ARMA model solutions developed in Chapter III follow 


19 





as a special case. Appendix I then presents two examples of 
nonlinear ARMA modeling. First a somewhat academic example 
of a cascade of linear and nonlinear subsystems is given 
then a nonlinear ARMA model is proposed for the tracking 
behavior of a phase locked loop. 

Finally, in Chapter V, two applications for the linear 
and nonlinear ARMA modeling methods developed in Chapters 
III and IV are discussed briefly. (They are reduced order 
modeling of complex systems and modeling for fault detection 
and diagnosis.) Then in Section V.B. conclusions are drawn 
on the results of this work and a list of significant open 
questions (both old unanswered questions and new ones 


raised here) is compiled. 
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ЖЕ  ISCRETE TIME LINEAR SYSTEM MODELING; BACKGROUND THEORY 


While few physical systems are absolutely linear, linear 
models often suffice to accurately describe their behavior 
under normal operating conditions, A rich body of theory 
has therefore been developed for the analysis and modeling 
of linear systems [Ref. 22] and a thorough knowledge of 
EE theory is vitally important to anyone interested in 
understanding the functioning of these systems. The con- 
tinuing expansion in the availability of powerful, inexpen- 
Sive digital computing capabilities has also made discrete 
time techniques take on a special prominence. With this 
mo tivation, the portion of the background theory in 
discrete time linear modeling upon which much of the re- 
mainder of this work depends, is developed here from the 
unifying standpoint of a minimum mean square equation error 
model solution. 

The moving average and the autoregressive models are 
developed first for single input single output systems. 
Their solution via the Levinson algorithm is presented and 
from this algorithm alternate solution methods based on 
Hattice filter structures are derived. It is shown that 
most all of these results can be generalized to the 
multiple input multiple output case and the corresponding 
multichannel modeling methods are developed. Finally, 


adaptive implementation of the modeling methods for both 


a 





the conventional filter structures and the lattice filter 
structures 1S presented as an alternative means of obtaining 


Medel solutions. 


A. MOVING AVERAGE MODELS 

The moving average (MA) model was among the earliest 
discrete models developed, [Refs. 4, 11 and 19] It estimates 
the current value of the output of a system as a welghted 
combination of the present value and N past values of the 
system input where N is the order of the model. The problem 
then is to estimate the weighting function or impulse re- 
sponse of the MA model in some fashion. Since the MA model 
characterizes a system in terms of a finite duration approx- 
imation of its impulse response and since any linear time 
invariant, single input single output ДЕТ 1s completely 
specified by its impulse response, the MA model is quite 
general and can be used for a wide class of systems. De- 
fining (N+1l)-vectors of model weights and input data as 


Т 


aa] > 


Аа) 





E ресет" г" 15 used to indicate that in spite of 
the fact that these vectors are used for a N-th order model, 
they are (N*l)-vectors with elements indexed from zero to N 
Eher than from one to N. Superscript T demotes transpose. 


Superscripts in parenthesis will later be added to the 
model coefficient vectors to explicitly indicate their de- 
pendence on the order of the problem being solved. They are 
omitted for simplicity however whenever doing so does not 
result in ambiguous notation. 
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ет 
ЖСК uk) = ulk-N)] (ООШ) 


the MA estimate of the system output becomes 
| Т 


у(х) = Y almdulen) = u (Kk) а 
n=0 


ae 


ems of the modeling approach of Figure 1.3, F and 


20 
Pag areassumed to be zero. Fi lS a linear time invariant 
function of past and present values of u(k). Assuming sta- 


tionarity, an expression for the mean square value of the 


E "asa quadratic function of the weights {а(п)} 13 


given by 
T ШІ 
+ + + 
Е, EU M cNP Lt та 2222 
u ul uy W 
where in general Roy) are), Em e ivCk)wCk*n)], 
ри. = e{v(k)w(k) +} and el } denotes expectation. 


В (0) 9: В. (-N) 
ча ша 
A hs 
U u 
в) 
uu аи 
Р (КС +--+ Rk сой 
o y uy 


pie surface described by equation 2.2 can be pictured as a 


concave hyperparaboloid with a unique minimum and the 


2 





characteristics of such a surface are described in Appendix 
Е. For example when N=1, the MSE as a function of a(0) and 


a(1) appears as shown in Figure 2.1. 


a(0)—» 


Figure 2.1. MSE as a function of model weights 
for a first order (N-1) MA model. 


The minimum mean square error solution for the coeffi- 


cients is given by 


R 
u 


+ +а орт ^ Z4 0223) 
u uy 
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(Пе corresponding minimum value of the cost function 


E ENSE) ca ig (2.4) 
У — 


Equation (2.3) is a discrete time matrix form of the Wiener 
Hopf equation 


со 


Bo (E f КОС OO (2.5) 
uy uu 


- СО 


meme u(t) and y(t) are the continuous input and output 
Signals and h(t) is the system impulse response. The pro- 
cess of finding Bm in equation (2.3)is the discrete time 
equivalent of  deconvolving the input autocorrelation func- 
Em irom the cross correlation of input and output to 
obtain the system impulse response in equation (2.5). 
Consequently the MA modeling process has been called dis- 
crete Wiener filtering or stochastic deconvolution. 

This model constitutes a direct form approach as defined 
ШО СО н Оп Т.1 but does not encounter difficulty in obtaining 
the model weights since both Fig and Pao are assumed to be 
Zero. As such, it possesses the advantage that the estimates 
of the model parameters will not be biased by the presence of 
additive noise on the output of the system as shown in 
Figure 2.2, as long as the noise is uncorrelated with the 
input signal. This can readily be seen by replacing y in 


equation 2.3 by y * v. Additive noise on the input signal 
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MA MODEL 


Figure 2.2. Moving Average Modeling 


however will adversely affect the modeling process, and 
introduce a bias in the solution for the model coefficients. 
In the transform domain, the model can be represented by 


a polynomial in powers of ae and has therefore been referred 


to as an all zero model 
N 
ICONE A a (2.6) 
met) 


Momcerms Of this transfer function relationship, any bias 
introduced in one or more of the model coefficients has the 
Beer Of Shifting the zero locations of the model. 
In summary a discussion of the advantages and disad- 
vantages of MA modeling is instructive. 
Advantages: 
1) The solution for the model parameters involves 


only linear equations. 
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The solution is unbiased in the presence of 
additive noise on the system output as long as 


the noise and system input are uncorrelated. 


3) Since the model is nonrecursive it is always 
сае 
Disadvantages: 
1) The number of terms (N+1) needed for sufficient 


2) 


ED 


4) 


model accuracy may be quite large. 

The solution of a large system of equations is 
required. 

The required correlation terms are usually not 
known and must be estimated by assuming ergodi- 
city and averaging in time. This requires the 
data to be windowed and set to zero outside the 
averaging interval in order to maintain the even 
symmetry of the autocorrelation functions. 

The modeling process is restricted to linear 


time invariant systems. 


Pe LOREGRESSIVE MODELS 


The autoregressive (AR) model attempts to deal with one 


the model. 


e(k) 


of the difficulties (1) encountered in MA modeling; the need 
for a large number of coefficients to accurately describe 
БЕ eel 19 and 28] In AR modeling, 
which is sometimes referred to as linear prediction, a pre- 


diction error approach is considered where 


N 
= y(k) - b(n) y(Ck-n) с Та) 
n=l 
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Brass can be written as 


a s WS S099 


= y(k) - y(k) 1b 2215) 
with 
T 
В OSI *** vCOK-N)] Gare) 
and 
T 
БИЕ ТЫС Боо] (Сета) 


Here Fio and Pao are assumed to be zero (this assumption 
Will be modified later to allow a dependence on the input 
Signal in the synthesis phase) and Poo provides an estimate 
of the current value of the system output as a weighted sum 
of N past outputs. The mean square value of prediction 
error as a quadratic function of the weights {b(n)} is 


given by 
r Terme) ова 
= y 


and the corresponding MMSE solution for the weights is 


given by 


A D = G2 8b) 
-УУ-ОРТ Y 
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E = (0) - e Ge 


а 
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Using equation 2.7a, an expression can be written 
in the transform domain for the prediction error model which 
accepts y(k) as its input and produces the error sequence 


ЕШ) ас its output. 


N 

EX) - -n 

YZY = l - > Im - B z) га) 
mel 


If it 1s assumed that the system input output relationship 


can be represented in transfer function form with 


иа). auo) E 
772) = H(z) ЕЕ = Biz) О) 
1- Dem) 2 
n=1 


and that the model parameters can be determined so that B(z) 
setz), then the prediction error output will be exactly 
e(k) = a(0) u(k). For this reason AR prediction error 
modeling has often been called inverse filtering since the 
prediction error filter essentially reverses the actions of 
the system (with the exception of a gain). Since the 
analysis model is in this inverse form rather than in the 
direct form, the presence of additive noise on the measure- 
Ment of the system output as shown in Figure 2.3 will intro- 


duce a bias in the solution for the model parameters. This 
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Figure 2.3. Autoregressive prediction error 
modeling as an inverse filtering 
process. 


has the effect of shifting the roots of B(z) which are es- 
timates of the poles of the system and is the price paid 
пе ability to obtain the model solution from a set of 
linear equations. 

Maas tar, only the analysis portion of the AR modeling 
process has been discussed. With the inverse filtering 
interpretation of the prediction error analysis model, a 
reasonable synthesis model is given in transfer function 


form as 


(CZ) = ВЕТ 





with the gain term set so that the mean square value of 
a(0) u(k) is the same as that of the prediction error 


Aenal. Thus it follows that 


ее: 


р (0% 
uu 


ao х= (2.12) 
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Since the synthesis model is in the form of an all pole 
filter, an appropriate impulse response with infinite dura- 
tion might be obtained using a low order model (small N), 
Ещиеви that is impossible to obtain in any finite order МА 
model. This is not to say however that a low or even finite 
order AR model will always be an appropriate model for any 
linear system, ТЕ the transfer function representation for 
E Siem contains both zeros and poles, no finite order AR 
SOMA model can serve to exactly represent it. This fact 
can be understood by considering the form of a geometric 


aries 


EE (2.13) 


о 
— m (cz 71) oa т 
1-Cz n=0 

which shows that a single pole can be represented Бу ап 

infinite number of zeros and visa versa. Thus if the sys- 

tem has a single zero, a high order AR model may be required 

ШІ Гергесзепі it with sufficient accuracy. 

In summary, the advantages and disadvantages of AR 
modeling may be listed as follows: 
Advantages: 

1) The solution for the model parameters involves 
only linear equations. 

2) Sometimes an appropriate infinite impulse re- 
sponse can be obtained with a small number of 
parameters in the model. 

3) Direct knowledge or measurement of the system 


input is not required for determing the system 


poles. Only a knowledge of its mean square 


a 





value is necessary for determining the gain 
ACEO 

Disadvantages: 

1) The model is biased by the presence of additive 
noise on the measured system output signal. 

2) The number of terms required for sufficient 
model accuracy may be quite large if zeros are 
бшсе ые system. if this occurs, the 
inversion of a large matrix will be required. 

3) The required correlation terms are usually not 
known and must be estimated by assuming ergodicity 
and using time averages. 

4) The modeling process is restricted to linear 
time invariant systems, 

This list of advantages and disadvantages is quite 
Similar to the one compiled for MA models with two notable 
differences; the bias in the model and the absence of a 
requirement for input measurements. This second point is 
significant in that it has led to the application of AR 
modeling to many problems where an input signal is unmea- 
surable or indeed does not exist including speech modeling 
and spectral estimation. [Refs. 2, 5, 12, 15, 21, 32 and 44] 
The noise problem has restricted the process to applications 
where measurements with sufficiently high signal to noise 
ratio are available, making the effects of the bias minimal. 


Meets. 24 and 43] 
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Seek eE CURSIVE IN ORDER SOLUTIONS FOR AR AND MA MODELS 

The preceeding disucssions of the AR and MA modeling 
problems tacitly assummed an apriori knowledge of the correct 
model order. If this knowledge is not available a reasonable 
approach for determining the correct model order must be 
developed [Ref. 53]. A commonly employed strategy is to 
successively increment the model order while observing the 
MSE until further increases fail to substantially reduce the 
MSE. This requires solving for a number of different models 
and can be an arduous task if equations (2.8b) or (2.3) are 
employed directly. 

The autocorrelation matrices appearing in the AR and 
MA model equations (2.8b) and (2.3) are highly structured 
matrices (both Toeplitz and symmetric) and this fact can be 
exploited to facilitate the solution of these equations. 
The Levinson algorithm (Refs. 9, 10 and 27] makes use of 
this structure to obtain model solutions recursively in 
order, that is, the solution for the n-th order model is 
assumed to be known and the solution for the (n*1)-st order 
model is then obtained from it. In this manner it is pos- 
Sible to start with a first order AR or a zeroth order MA 
solution given by a single equation and build up the 
desired order solution. The AR model will be treated first 
since it is a special case that simplifies the analysis. 
The simplifications arise due to the fact that the Es vector 
on the right hand side of equation (2.8b) is made up pri- 


marily of terms also appearing on the left hand side in — 
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Superscripts in parenthesis are used to explicitly indicate 
the order of the problem when specifically needed. 
1. The Levinson Algorithm For AR Modeling 

The n-th order AR model solution of equation 2.8b 

is given by 
а м (п) (2218) 

ШУ му у 

The Levinson algorithm assumes a relationship between the 


n-th and (ntl)-st order solutions given by 


(та) КО 


D 
E = э i E (2.15) 
Q x, (n*1) 
and solves for the vector un and the coefficient ky nt, 
Define permuted versions of the vectors Ба) апа Es as 
ща) апа ше by reversing the order of their elements. 
1 
ben?” а 
ЕП = o | Ge) 
- E ; =y y 
n 
yy | 


Because of the Toeplitz symmetric structure of the auto- 


correlation matrix, equation (2.15) can also be written as 


(ба) (ЕЛ | (п) 
Avy Б = By ЕШ) 
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and this relationship is essential in the development of the 


Mevinson algorithm. 


(To apply the algorithm therefore when 


time averaged estimates of the needed correlations are used, 


the data must be windowed prior to averaging to maintain the 


even symmetry in the autocorrelation function estimates and 


produce the required structure in 
Making use of equation (2.15), in 


of equation (2.14), and using the 


solve for eu and ntl) results 
m _ г 
and 
Ш 
(п) 
k = 
ROOTED 
yy D уу 


ЕЕ, in using equations 2,15 and 2.18 to obtain b 


en 


From b 


information, „Intl 


, need be calculated. 


the autocorrelation matrix. 


the (n*tl)-st order version 


polae onship of (2177 to 


aol 


ОЩЕ ga) 


(n) 


— 


ЕВ) 


en) 


(n+1) 


via the Levinson algorithm only one new piece of 


The denominator 


of equation (2.18b) can also be recognized from equation 


(2,8с) as the MMSE for the n-th order AR model E 


mS there is little concern over 
ШО, Lf the n-th order solution 


tion (zero MSE) there is no point 


prediction by increasing the order to ntl. 


Ga 
2 9 


CMe DOSSIDILIEY GE it. being 


and 


produces a perfect predic- 
ес отит ЕО на га better 


The evaluation of 


equation (2.18b) can be further simplified by observing that 


the MSE also follows a recursion from one order to the next 
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given by 


о _ Ge 


2 = E, [1 - К 


E ] 290 
making it unnecessary to evaluate the denominator at each 
Eus of n. (Details of this derivation are omitted here 
but included in Appendix A in the derivation of the more 
general multichannel Levinson algorithm.) This relation ino te 
the propagation of mean square prediction error also leads 
Бе ће restriction that a inti) must be bounded in magnitude 
ШЫ unity. 

2. The Levinson Algorithm For MA Modeling 


Next consider the n-th order MA model given by 


(п) 

ca) + 2 (n) 
-—— a en 02220) 
u uy 
and again, assume a relationship between the (n*1)-st order 


and n-th order solutions given by 


(n) 


и ml) 
„(n+1) 2 H 
a = ¡(e + E (22010) 
0 Mc 


— 


Meerce that in this n-th order problem, R E 


aeu 
(ne) 
И : 


is actually a 


Бу ntl matrix and could be written as К, Define 


36 





- 
сие eee PR, 61) | CHE Em 


R (пт) Dane 
uu 


Using equations (2,21) and (2,22) in the n+tl order MA model 


equation it follows that 


(n+1) (ntl) (ntl) 
Y zu g 


сар 


ов а) 


where nr) 


is defined in a manner similar to (2.16) and is 
Шипргісесй о: the coefficients that arise in an (ntl)-st 
order autoregression on the input signal u(k). 

Therefore to obtain a moving average model relating the 


System output signal y(k) to the input signal u(k), an 


autoregressive model for the system input must first be 


solved. Furthermore, 
Dm) 
+ 
eu К, шы er а 
g = MM (2.23b) 
в (0)- (ТЕ) gintl) 
uu цц — 


ea the denominator of equation (2.23b) is the MMSE in the 
(ntl)-st order autoregression on the signal u(k). 

MS зара в теате to note that in applying the 
Levinson algorithm to find a given order AR or MA model, all 
lower order models along with their MSE's are obtained. 


(п) 


Also, intermediate quantities emerge (the{k } in the AR 
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model and the П апа T 


) in the MA model) which fully 
eharacterize the models and could be used as an alternative 
Berne ia(n)} or tb(n)} coefficients. This point will be 


developed further in subsequent sections. 


NA TTICE FORM AR AND MA MODELS 

The Levinson algorithm derived in the previous section 
can be used to derive lattice structures to implement the 
MA model and the AR analysis and synthesis models as alter- 
natives to a tapped delay line type of implementation using 
BE Ncoefficients a(n) or b(n) directly. [Refs. 29, 30, 32 
end 33] 

1. The AR Modeling Lattice Structures 

From the relationship between the (n+tl)-st and n-th 

order solutions to the AR modeling problem determined in 
Sama tions (2.15) and (2.18a) it follows that the transfer 
function of the prediction error model can be written re- 
cursively in order as 


BEDI) e2) я ВО) , y (nt D,-n-1g() (71, (2.24) 


Defining a new transfer function 


Bz) OS (2298) 
equation (2.24) can be written as 
dd ue int ign (2.255) 
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and an expression can also be written for Br recur- 
sively in order by rewriting equation (2.25а) for order 


О апа substituting equation (2,24) yielding. 
دا 0( )0 د‎ ۹ орк В? Вет) (а) (2.250) 


Pemoiscussed earlier in connection with equation (2.9), 
ВП? (=) describes the n-th order prediction error model and 
СП 155 input is the system output Y(z), it produces the 
n-th order prediction error signal. 


gy Е gn) 


Cz 0 eJ (226) 
In the time domain this signal can be interpreted as the 
error in predicting y(k) forward in time from a 

weighted combination of the n past values (y(k-1)':''y(k-n)). 


(п) 


To understand the significance of B (z) consider the out- 


put signal when this model is excited by Y(z). 


gy = CIC) 


1) 


Ж 0 s P oz wo (2.27) 
L=1 
In the time domain, u can be interpreted as the error 


in predicting y(k-n) backward in time from a weighted 
combination of the future signals fy(k-n+l):-‘‘y(k)}. These 


n-ath order forward and backward prediction processes at time 
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Ko E iiiustrated in Figure 2.4. Henceforth, an overbar 
Will always be used to denote quantaties associated with 


backward in time predictions. 


yao 4 
| | 
| Forward | 
| Prediction | 
[ Е SSS | 
| Backward | 
Prediction 
--| ы... | 
| | 
| 
| 
| 
| 
| | 
B - 
Kg-n Kg k 
Figure 2.4. Forward and Backward Prediction Error 
Filtering. 


ШЕП equations (2.25b) and (2.25c) equations can be written 
recursively in order for these forward and backward predic- 


tion error sequences аз: 


Mu!) nb _ 31) (п) 


(К) = е (К) (k-1) O77 28a) 


E (n) Е „(n+1) (n) 


Е (Rei) (k) 265) 





Noting that the prediction error for a zero-order AR pre- 
E Sr of yk) Cor no predictor at all) is just the signal 


Eu itself, 


Din (2.280) 


jen prediction error filter can be drawn in lattice form as 


shown in Figure 2.5 for a second order model, 





51) (р) 502) (р) 


Figure 2.5. Lattice Form ОҒ A Second Order Pre- 
emeti on Error Model. 

This structure has many interesting properties, 
among the most important of which is the successive de- 
ling property. In going from one order AR model to the 
next, all of the previously determined transfer function 
coefficients {b(n)} will generally change. The Levinson 
algorithm shows however that only one new piece of infor- 
mation is needed to obtain the optimum (n+1)-st order 
solution from the optimum neth order solution (see equation 


Ш 21)). Іп terms of the lattice filter of Figure 2.4 this 





means that given the optimum n-th order model in lattice 


morn, one need only add another stage to the structure, 


(n+1) 


setting the coefficient of that stage k Воо се 


the mean Square value of ао 


(к). Nothing in the first 
п stages need be changed, The overall high order minimi- 
zation problem is in this fashion decomposed into a se- 
quence of first order minimizations, one at each lattice 
сасе. 

Another important property of the lattice which can 
be proven and will be of use later is the orthonogalization 
of the backward prediction error sequence [Ref. 32] which 
EHues that 

0 123 
-(3) 


ШЕ (oe Cas) = 02509) 


(i) M" 
2 Eus 


tij 


Thus it is seen that a set of orthogonal signals (the back- 
ward prediction errors at the various stages) are generated 
as a by-product of the lattice model. 

As a consequence of the successive decoupling pro- 
perty of the lattice, a number of alternatives to equation 
(2.18b) for determining the lattice coefficients can be 
found. The most obvious method is to set Karl) to 
explicitly minimize the mean square value of forward pre- 
diction error in equation (2,28a) at the (n+tl)-st order 


stage given the best lattice of order n, This is termed 


the forward method and is denoted by a subscript F on the 
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ШЕГІСе coefficients. The resulting solution is given by 


Greta SUR) cy 153 


k - 025.515 
ЕА 


Р 
Alternately, the mean square value of the backward predic- 
tion error signal in equation (2.28b) could be minimized to 
determine the coefficient resulting in the backward method 


solution given by 


(n*1) MO US x 


k = C23) 
(л) уу 27 


Р se 


Since, however, the forward and backward prediction error 


transfer functions are given by g 0) (25 and RUM 
it follows that 
в (еур = [BO (| (2.32) 


and since they are both driven by the same input, Y(z), the 
mean square values of both the forward and backward predic- 
tion error signals at a given stage are the same making 
equations (2.30) and (2.31) equivalent. It is also possible 
to show that they are equivalent to equation (2.18b). 
Recognizing that the required expectations will eventually 
have to be estimated by using time averages, these two 


methods for calculating ae 


will not in general be 
exactly equivalent and it might be preferable to use the 


arithmetic mean of the mean square values of forward and 





Beeiward prediction error as a cost function 


ЕЕ A y (2.33a) 


This leads to a third method derived by Burg [Ref. 5] in 


his work on maximum entropy spectral analysis and given by 


(ntl) _ а G1) } 


в 


Пе) 7 
BG snc) 
NO се that EU is the harmonic mean of pe and pus, 


ENS Pth method due to Itakaura and Saito [Ref. 20] can also 


ШЕ. бегіуес which results in 


(п) 


6:9 = 
15 т СС 2 
€ fe e cs Ml e 
(n+1) е « е 
апа Krs is simply the geometric mean of the forward and 


backward coefficients. 
Since equation (2.34) is of the form of a normalized 


correlation k will always be bounded by unity in magnitude 


iS 
as required by equation (2.19). Furthermore since 
| Harmonic Mean | < | Geometric Mean | 
Gai) 


it follows that k will be similarly bounded. These 


BG 
bounds are significant since Markel [Ref. 32] has shown that 


(п) 


Ik" "|«1 is a necessary and sufficient condition to ensure 





Ec the roots of gin) 


(z) be within the unit circle guaran- 
teeing the stability of the n-th order all pole model. No 
such guarantees of model stability exist when the forward 
or backward solution methods of equations (2.30) and (2.31) 
are used with the correlation estimates obtained by aver- 
aging for finite time intervals. 


To determine the AR synthesis model in lattice form 


it is only necessary to rewrite equation (2.28a) as 


(n) Саа) ОПЕ) —(п) 
е е е 


(k) БК (kel) 02.35) 
Together with equations (2,28b) and (2.28c), this describes 
Mies tructure shown in Figure 2.6 for a second order case 


and when it is driven by the second order prediction error 


Seema), it will reconstruct y(k) exactly. Thus it imple- 


ments the transfer function 5 or, in general, when 
B CZ) 
1 Xm : 
stages are used --7= -. . Recognizing that if the pre- . 
ВС) (а) 


diction error model is an accurate model of the system 

denominator polynomial, И». ium De sme 
is used in the synthesis model. Because of analogies with 
transmission lines and wave propagation models, the lattice 


Mu scjients have been referred to as reflection coefficients. 
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Бігіге 2.6. Lattice Form Of The All Pole Synthesis 
Model For The Second Order Case. 


EE A Modeling Lattice Structure 
A similar lattice form is applicable to the MA 
modeling problem. From equations (2.21) and (2.23a) the 
transfer function of the MA model can be written recursively 


order as 


ent) (n) (ntl)sn+l 


А (2) = А`П’(2) + g (2) 9050 


where, as discussed in connection with equation (2.23а), 
m 


B (z) is the backward prediction error transfer function 
for an autoregressive model of the input signal u(k). 
Pultipiying both sides of equation (2.36) by U(z) and trans- 
mine into the time domain it follows that 


ет 1 


O y o + р 1760р) фа зи 
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+ eee E 
E Lk) ШЕЕ раната prediction error signal from 


where е 
the autoregression on the input signal u(k), and can be 
Obtained by operating a prediction error lattice with u(k) 
Ве input. Then with the additional term in equation 


EN the lattice form of the MA model can be drawn as 


Ве їп Figure 2.7. 
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г: 
E n~ EO) s P^ qo 
(0) 
Е „m RS 
y 9^ do y ر‎ y ^^ do 


 епге 2,7. Lattice Form Of The MA Model (Second 
Order Case). 

It was stated earlier that the AR prediction error 
lattice, as a by-product, forms a set of orthogonal or 
uncorrelated backward prediction error signals from its 
input. Here in the MA model, these orthogonal signals are 
linearly combined to form the MA estimate of the system 


But. If the input signal u(k) is a white process, ап 
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examination of any of the solution methods previously dis- 


(п) 


ВЕ will show that all of the {k PS Е ав азс ficients 


will be zero since the delayed samples of a while process 
—ы—— — 


(xj | 


E ready orthogonal. Otherwise the {К еее 


coefficients will be set to orthogonalize the backward pre- 


E ON error signals. As a consequence of this, the 


Grp 


weighting coefficients g can be set independently of 


бы) 


each other; that is g can be set to minimize 


p 


М ЕЮ уп” (2.38) 


en) 


given the best prediction of order n-l, y (к). This results 


(n) 


in an alternative expression to equation (2.23b) for g 


given by 
па —(п) 2 
„m _ le "OO 97700) — yoo E o) "m 
eM GO al) 
Here еп? (к) is the error between the system output and its 


n-th order MA estimate. 


Е, MULTICHANNEL AR AND MA MODELING 

Both the AR and MA modeling problems previously dis- 
cussed, as well as their solution via the Levinson recursion 
and lattice filter methods, can be generalized to the 
multichannel case by replacing the various signals with 
signal vectors and replacing the weighting coefficients 


with appropriately dimensioned matrices of coefficients. 
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ИС С155іог оГ this appears in Robinson [Ref. 45]. The 
Tons that describe the AR and MA models and thelr MMSE 


solutions are repeated here for convenience. 





Р № 
АВ о в b(t) y(k-i) (2.40a) 
Tal 
ЕКИ (о) о ba) в © 
ЖУ ey уу УУ 
ЕШ (1) 54092220201 Қ ШӨЛІ 152) БІТІСІ) 
yy уу Yy = yy 
2 : (72405) 
В в (N-2) .... R. (0) b(N) к (М) 
К 502) уу | ЕШ 
^ N 
MA NP a) ue) (2.ч1а) 
1=0 
R60) касы. REN) a(0) RC? 
D) Ru R 51-00 a(l) И В 
(2,41b) 
В (М) КОМ ООШ ЖЕ LL в (0) a(N) в (М) 
uu uu uu ШЕЛ 


In a multichannel generalization, y(k) becomes a vector of 
00 output signals, u(k) becomes a vector of 0. input signals, 
b(i) becomes a square matrix of % х 9, coefficients and a(i) 


becomes a Qo x 0. Matrix of coefficients so that the N-th 





order multichannel model equations can be written as 


А N 
AR у(к) 22, b(n) y(k-n) (2.42) 


n=l 


MA y (k) a(n) u(k-n) (2.425) 


0 


a Mz 


EN NScuations for the MMSE solutions (2.40b) and (2.41b) 
generalize directly as well by replacing each correlation 
coefficient ос) by matrices of correlation coefficients 


given by 
Bann) - Е СЮ) Ол) 02:13) 


ЕВЕ Vik) and w(k) are signal vectors. This causes the 
overall correlation matrices to take on a block Toeplitz 
structure. The transfer function relationships of the AR 
prediction error model and the MA model take the form of 
matrix polynomials 


il 


Ьб1) 271 + +... + Ъ(М) 27 


В(2) Gala 


1 


А(2) = а(0) + а(1) 27 + +++ + a(N) 2S 


(2.44b) 


БО that 


E(z) = EL - B(z)] т) (2.440) 





Y(z) < A(z) U(z) (2. чна) 


where Е(2) is the transform of the multichannel AR prediction 
1 vector e(k) = y(k) - y(x). Alternately, equations 
(2.44a) and (2.44b) can be written as single matrices whose 
ES are polynomials in Z rather than scalars. B(z) is 

of necessity a Oo p 09 Square matrix polynomial while the 
dimensions of A(z) (Qo x 0.) depend upon the number of 

inputs and outputs which need not be the same. 

In the single channel AR problem, B(z) provides the 
transfer function of the prediction error, or inverse fil- 
ter, and must be inverted to obtain the all pole synthesis 
filter. The stability of the synthesis model therefore 
depends upon the roots of this polynomial. The matrix 
polynomial [I-B(z)] in the multichannel AR problem is, in 
like fashion, an inverse filter representation and must be 
inverted to obtain the synthesis model. This inversion of 
a matrix with polynomial entries is defined in the same 
manner as the inversion of a square matrix with scaler 
entries, To see what this inverse matrix polynomial looks 
like consider as an example a two channel autoregression 


With a prediction error filter given by 


5 cz) =b, ,(Z) 


TI 1:2 


ЕЕ 2) В ер Са) 


SI 





Applying Cramer's rule, the inverse matrix polynomial is 


written as 


1-5,,(2) 5,762) 
-1 _ ald 
ГІ-В(2)) = det[I-B(z)] 
5152) 1-51 (2) 


НЕСІ 1С is apparent that the stability of the multichannel 
synthesis model is dependent upon the locations of the zeros 
ШЕ (Пе polynomial det[I-B(z)]. 

This straightforward generalization of the AR and MA 
models is what has customarily been used in the literature 
to develop the multichannel models and similarly derived 
generalizations of the Levinson algorithm to solve them 
recursively in order are available as well. [Refs. 57 and 61] 
The multichannel AR and MA modeling problems and their 
solutions via the Levinson algorithm can however be recast 
as shown in Appendix A to make them compatible with the form 
of other linear and nonlinear modeling problems. To avoid 
confusion later in the application of the results, the 
derivations in Appendix A have been carried out in a generic 
form with x and d used to represent some of the signals and 
coefficients respectively. The symbols u, y, a and b have 
been reserved to denote system input, output and weighting 
coefficients. 

Baeetons €A.7) and (A.26) provide the MMSE solutions to 


the multichannel AR and MA modeling problems in forms 
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different than (although entirely equivalent to) those 
resulting from the straightforward generalizations of 
equations (2,40b) and (2,41b), The multichannel generali- 
Sn Of the Levinson algorithm derived in Appendix A 
can, with one exception, be seen as a matrix algebra 
generalization of the single channel algorithm and as, one 
would suspect, the single channel algorithm results as a 
Special case of Appendix A. The one exception is that in 
multichannel case, the n-th order forward and backward 
prediction error models are not simply related to one 
another. The single channel AR backward predictor is given 


-пвбп) (2-1; but in the multichannel case the backward 


use 


Dy Z 
Meedietion is not аар e 
: Because of this, two reflection coefficient matrices 
K and K are required at each stage in the recursion to 
relate the n-th and (ntl)-st order solutions rather than 
just one as in the single channel problem. Also, in the 
single channel case, the fact that ja? (а) | < 18722) | 


a dE made possible the de- 


and therefore efe 
rivation of the Burg algorithm and the Itakaura-Saito al- 
gorithm, which ensured the magnitude of the reflection 
coefficient was bounded by unity and that equation (2.19) 
Boma result in nonnegative values of MSE. In the multi- 
channel algorithm however, the forward and backward predic- 
MOL error covariance matrices and their traces are not the 


(0) апа в did as a result, straightfor- 


Ш ІШЕ (except for P 
ward generalizations cf the Burg and Itakaura-Saito algorithms 
are not possible. 
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Consequently with correlation estimates obtained by averaging 
over finite time intervals, there are no guarantees that 
equations (A.18a) and (A.19b) maintain the positive definite- 
ness of the prediction error covariance matrices. 
Multichannel generalizations of Burg's algorithm due to 
wei [Refs. чо апа 511], Morf [Ref. 39], and Strand 
[Ref. 50] which guarantee the positivity of the covariance 
matrices. are available but are not explored here because. 
of their complexity and because they would take the dis- 
cussion too far afield. 
The relations of equations (A.20) and (A.30) which are 
repeated here for convenience permit the construction of 
the multichannel AR analysis and synthesis lattice 
structures and the MA lattice structure. For the multi- 


channel AR model, 


Jt 
E o O (2.453) 
-(n+1) —(п) -(n+1)" _(n) 
Өрт у зш оек al Ge (2.45b) 
ЕО -705)-х0 (2.450) 


ШІСІ the corresponding prediction error lattice is shown in 


ENSure 2.8. 
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БО! (ко а’. 

Figure 2.8. Multichannel AR prediction error lattice 
SEI UCU [OF ada second order model. All 
signal paths are vector paths and summa- 
tions are vector summations. The multi- 


ра вота Indicate premiltiplication of 
the signal vector by the specified coef- 
iM етеп тарах. 


To obatin the multichannel AR synthesis lattice, equation 


(2.45a) can simply be rewritten as 


DE ай ) 


So + K 


(k-1) (2.46) 


ШОШ у пес in the structure shown in Figure 2.9. 
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|. Figure 2.9. Multichannel AR synthesis lattice 
structure for a second order model. 
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For the multichannel MA model, equations (2.45) and 


m) 


E A 
= y pau on 05 (2.47) 


ogg (к) +6 


describe the lattice structure shown in Figure 2.10. 


С 


Figure 2.10. Multichannel MA lattice structure 
for a second order model. 


As in the single channel case, the multichannel predic- 
tion error lattice exhibits the successive decoupling pro- 


perty and orthogonalizes the backward prediction errors at 


the various stages so that 
| | 177 
aa). ol: (2.48) 


{СО 


ШУ] 
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As a consequence of the successive decoupling, the forward 
and backward reflection coefficient matrices at the (n*1)-st 
stage can be set to minimize the trace of the forward and 
backward prediction error covariance matrices respectively, 
given the best lattice of order n. This provides alterna- 
tives to equations (A.17) for calculating K and K and also 
generalizes the forward and backward single channel solutions 


discussed previously, resulting in 


zu T 
oe и Б ща (2 цоа) 
о (191 (п) 
K =P A VETIS 
where 
o^ uere aget? ces T (2.1490) 


It is also possible to show that these relationships are 
entirely equivalent to equation (A.17), In the multi- 
channel MA lattice, the orthogonality of the backward pre- 
Eon еггог signals also allows the G matrices to Бе 
calculated in succession providing an alternative to equation 


(A.28b) and generalizing the single channel solution given 
(п) 


Myeequation (2.39). Setting G to minimize the trace of 


fe error covariance matrix 


Ga), (n) 


(a) 
E 0 


Ро = ele, Ce T ER 022502) 
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) where 


Sem) 


дут (2.50) 


е 
—0 
results in a solution given by 


im E 


G Е ра 


є (=) (к) у(к)1} (2.50c) 


Another important characteristic of the lattice solutions 
to the AR and MA modeling problems given by equations (2.59) 
and (2.50) and their single channel counterparts is that they 
do not impose any requirements to window the data when finite 
time averages are used to estimate correlations. The auto- 
correlation function of a signal is inherently an even func- 
tion so that и) = Bar This fact is responsible for 
much of the special structure of the correlation matrices 
that appear in the model solution equations, and was also 
used in the derivation of the Levinson algorithm, In esti- 
mating the autocorrelation function via time averaging over 
a finite interval, a window function that is nonzero over 
only a given interval must be applied to the data to retain 
this even symmetry property in the estimate. If the data is 
not set to zero outside a given interval, end effects will 
destroy the symmetry so that К Сп) Ё R En) as depicted 


ИШЕ Т риге 2.11. 








Averaging Interval | 


Figure 2.11. Time averaging to estimate 
CE NOUE windowing. 


he lattice solutions of equation (2.49) 
is no requirement to make such an artificial 
about the data (that it is zero outside some 

These properties of the lattice solution 


Pesponsible for their initial use by Burg in 


v(t*n) 


correlations 


however, there 
assumption 
interval). 
methods were 


his work on 


maximum entropy spectral analysis [Ref. 5] and have con- 


tinued to generate interest in the application of lattice 


methods to other types of modeling problems. 
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Е. ADAPTIVE MODELING 

The LMS adaptive algorithm provides a well known alter- 
native method for obtaining the solution to the AR or MA 
modeling problems which does not require the estimation of 
Memmelations or the inversion of a matrix [Refs. 58, 59 
and 60]. This algorithm updates an estimate of the model 
solution vector at each time instant by an amount propor- 
tional to the negative of the instantaneous gradient of the 


Meee function; i.e., in a MA model, 
+ + 
a (k*1) з а (k) =- p Y (k) (2.51a) 


where u is a proportionality constant or adaptive gain. 
Since the actual gradient is usually not known, it is 
approximated by using the square of a single sample of the 


ror as an estimate of the MSE so that 





А 2 
V OO т E = -2 MES e(k) (222515) 
i а са (к) 
апа 
a^ (k*1) = a (k)+2u ч (К) e(k) те) 


In each of the models considered here, the cost function 
Шот trace P) is a quadratic function of the model weights 
and defines a concave hyperparaboled with a unique minimum. 


Ihe functioning of the LMS algorithm under these conditions 
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can easily be understood by considering the scalar case of 
ion (2,51) illustrated in Figure 2,12. 

As this illustration shows, the algorithm can actually 
diverge for too large a value of adaptive gain. The rate 
of convergence is also dependent on the size of the adaptive 
gain. Widrow has shown that for stability, the gain must be 
set so that 


к< ү < (2. 52а) 





max 


While, in the mean, the weight vector converges with an 


exponential time constant of 


ee 


еа 


П CSB) 
where A_. апа л are the smallest and largest elgenvalues 
min max 


Ne input autocorrelation matrix R , ,. From the stand- 


U U 
point of stability, U should be made relatively small but 


for rapid convergence, equation (2.52b) dictates that it 


should be large. Setting 
u ку (2.53а) 


where à is a normalized gain and 0 <a <1, equation (2.52b) 
becomes 


А 
za c max E 
T = E E ЮЗЬ) 
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dont ee 


a) Small adaptive gain; steady convergence toward the solution. 





b) Intermediate adaptive gain; oscilatory convergence toward 
the solution 


MSE 





c) Large adaptive gain; divergence away from the solution 


Figure 2.12. Behavior of the LMS adaptive algorithm 
for various adaptive gains. 


and for a wide disparity between the largest and smallest 
eigenvalues Ce << а’. convergence will be quite slow. 
This consideration becomes increasingly important when high 
order model solutions are obtained adaptively since the 
dimensonality of the input autocorrelation matrix will be 
high and the possibility of a wide eigenvalue disparity 
шиеатег. 

These same adaptive techniques have been applied by 
ит ънс [Refs. 16, 17, 18 and 311, to the AR and MA 
lattice filter structures derived in the last section. The 
key difference between the conventional adaptive filter and 
the adaptive lattice is that in the lattice, the adaptation 
ls carried out on a stage by stage basis for each of the 
reflection coefficient matrices while in the more conven-' 
tional approach, the entire weight vector is adapted, It 
has already been established that the lattice structure 
makes the model solutions recursive in order. Implementing 
the lattice adaptivity makes the solution recursive in time 
as well since the estimate of the solution at each instant 
is dependent upon prior estimates of the solution, 

The conventional adaptive filter algorithm forms an 
error signal as the difference between scme desired signal 


and its estimate; i.e, 


e(k) = у(к) - а" (к) Ти ck) (2.54) 





where y(k) is the desired signal, а“ Dc c ТЕМЕ vector 
and и" (к) ое тестот, апа the time update for the 
weight vector is given by equation (2.5lc). To derive the 
adaptive AR lattice consider equations (2.45) for a single 
stage. The lattice in general has vector error and desired 
Signals and coefficient matrices as opposed to scalar error 
and desired signals and a coefficient vector in equation 
(2,54) but such a generalization is straightforward. Com- 


Meeame equation (2.45a) to (2.54) it is clear that: 


Ш) Ао is analogous to the error signal; 
2) ао is analogous to the desired signal; 
3) Sk) 1s analogous to the input signal vector. 


(ntl) as a cost function and applying 


Meee the trace of P 
a LMS adaptive algorithm to determine the forward reflection 


Meerricient matrix it follows that 


(тк) 


ao 


gum DASS В 


к(п+1) (п )—(п) 
= = 


(k) + 2u (k-1)e 


7. 55а) 


With these analogies, equations (2.51c) and (2.55a) are seen 
to be virtually identical with the exception that (2,510) 
uses a scalar error to adapt a weight vector and (2.55a) uses 
a vector error signal to adapt a coefficient matrix. 

Proceeding ina similar fashion with equation 2.45b it is 
Clear that: 


Eris 


1) (k) is analogous to the error signal; 


2) E is analogous to the desired signal; 





B) a ото che input signal vector. 


. т, T : Т 
DEN the trace of ae 1) Еее FUNGCELON, the time 
update relation for the backward reflection coefficient 


MATIX iS 


g(n+l) (n*1) ED (n) 


(к+1) = К (k) + 20 ON qm 


550) 
For a MA lattice, equation (2.47) must also be considered. 


Multiplying both sides of (2.47) by minus one and adding 


IM results in 


Т 
ал ск) = en 7? (к) E з; (2.56) 
Gre! > ; 
where eg (k) is defined as in equation (2.50). It is 
evident that; 
co) | 
в) ез (k) is analogous to the error signal; 
2) е () 1s analogous to the desired signal; 
3.) Sud is analogous to the input signal vector, 
With the trace of E as a cost function, the time up- 
Gas) 


Ee relation for G is given by 


Бе) (131) СК) A ЕТ SS (n*1) T 


(к+1) = С (к) 


.-. 


(MEE 


257) 
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Ше 15 сіспі ісатт го note that three different adaptive 
gains have been used in equations (2.55) and (2.57) and 
that the gains have been superscripted indicating that they 
vary from one lattice stage to the next. For stability consi- 
derations the adaptive gain used in the LMS algorithm must 
satisfy equation (2.52a) and therefore is related to the 
largest eigenvalue of the input autocorrelation matrix by 
equation (2.53a). In developing the time update relations 
for the lattice coefficients, three different input signals 
were used and these inputs also differ from one lattice 
stage to the next. Indeed, even for the case where the іп- 
КОО О К) to the lattice structure is stationary, the inputs 
to all lattice stages except the first are nonstationary 
Since these inputs are outputs of other lattice stages. This 
fact indicates that time varying adaptive gains are appro- 
priate as well in equations (2.55) and (2.57). Equation 
(2.53a) is of no direct usefulness however in setting the 
adaptive gains since the time varying eigenvalues are not 
known. Recognizing that the largest eigenvalue is always 
less than the trace of the input autocorrelation matrix 
(which is a measure of the power in the input signal vector) 
the gains can be set as 


ul (x) S (2.58а) 
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=— > seb) 
‘+1 (К) 
(at ID) ek) = = messe) 
ы ые 
та: 
where o cor с Ue and y E ape estimates of the 
mut : me 1 ml 


power in the three input signal vectors and a is a normalized 
 Гіуе step size with 0 < с <1. A method that has commonly 
been applied to obtain these estimates is to employ a first 


order low pass filter so that 


(m) Т-(п) 


1+1 (+1) ° * E1-aJo „Со Ca (Kl) e (к-1) (2.598) 


(п) 


(к) 


ыл; (2.59Ь) 


к 2. - 2 
ШОО tl) = [1-alo,,, Ck) *ae 


(n) 


Y. 40024 2 E1-ody, LOO ^ fae, I? G0 Te СО (o (2.590) 
MEN together, equations (2.55), (2.57), (2.58) and (2.59) 
define the adaptive solutions for the AR and MA multichannel 
lattice models. 

To understand the potential advantage offered by the 
adaptive lattice form, recognize that while the conventional 
approach solves a high order minimization problem by adapting 
all the coefficients at once, the lattice breaks the problem 


Mew Into a succession of lower order minimizations at each 


stage and solves these lower order problems adaptively. The 


ст 





MemenslOnality of the input autocorrelation matrices. at 
each lattice stage in general is significantly less than 
chat of the large input autocorrelation matrix in the con- 
ventional adaptive algorithm and consequently it is hoped 
that the possibility of a large eigenvalue disparity with 
its attendant slow convergence is reduced. 

This advantage is most evident in a single channel 
adaptive lattice where the inputs at each stage are single 
Signals and their corresponding autocorrelation matrices 
are 1x1 in dimension. In this case the ratio of smallest 
to largest eigenvalues is unity and the convergence of each 
stage 1S quite rapid while the convergence of the overall 
model is independent of the eigenvalue ratio for the over- 
all higher dimension input autocorrelation matrix. This has 
been demonstrated by Satorius [Refs. 46 and 47] who has 
shown that the single channel adaptive lattice converges 
much more rapidly than the corresponding conventional adap- 
tive filter, and does so independently of the eigenvalue 
вето On the overall input channel autocorrelation matrix. 

Furthermore in a single channel adaptive lattice, the 
time update relations are simplified by the fact that the 
forward and backward reflection coefficients are the same. 
Using the average of the mean square values of backward 
ema rOrward prediction errors as a cost function and applying 


an adaptive algorithm it follows that 


AR 





ШІ) yeu. A | > 


С 1 СК) 


Cae) 


TS (k) J (2.608) 


where 


m) 


2. Оа) Ze 2 
9 ui (kt) -L1-a]o, , CK) +5Le (k) *e K (2.60) 


The very nature of the lattice structure however with 
the output of one stage providing an input to the next stage, 
greatly complicates the analysis of the convergence proper- 
ties of the adaptive lattice model. Even when x(k), the 
mpc со the lattice, is stationary, inputs to all stages 
except the first are nonstationary. An approximate analysis 
of convergence and stability on a stage by stage basis is 
possible if it is assumed that all prior stages have con- 
verged and are providing stationary inputs to the stage 
under investigation. With this assumption, the adaptive 


(ntl) E — anal) 


millon for the К matrices аге 


obtained from the operation of three independent LMS algo- 


—(n) 


rithms as shown earlier, with inputs given by e (к-1), 


—(nt T DES 
ДА ск) апа an mn еи ее limits on 


the adaptive gains used in the stage and the convergence 
properties of the stage are then determined by the eigen- 


p(n, drum E puni 


lies of the Р matrices. A more exact 


analysis of the properties of the adaptive lattice that 









considers the nonstationary character of the inputs to the 


second and subsequent stages is not currently available. 
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A ARA MODELING 


One of the most serious disadvantages of either AR or 
MA modeling is the fact that to adequately represent even 
Simple linear systems, both methods may require a large 
number of parameters (a high order model). This problem 
arises since, from a transfer function standpoint, AR and 
MA models attempt to model the system using only poles or 
only zeros, in spite of the fact that the physical system 
may have both zeros and poles. While modeling the effects 
of a zero with a number of poles and visa versa can be 
analytically justified as shown in the previous chapter, 
it makes far more sense (both from the viewpoint of model 
accuracy and efficient use of model parameters) to let the 
model represent the system as it really 1S with both zeros 
and poles if this is at all possible. The ARMA model is a 
generalization of the AR and MA models and accomplishes 
exactly this, representing the system in rational transfer 
Emetion form. 

MaS worth noting that the titles of all pole and all 
zero modeling that have been associated with AR and MA 
modeling are misnomers. Both haye equal numbers of 
zeros and poles. In the AR model however, all the zeros 
occur at the origin of the z-plane as do the poles of a MA 


model. The ARMA model removes these constraints. 
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After a brief discussion of two alternate ARMA modeling 
methods due to Shanks and Prony, the equation error formu- 
lation for ARMA modeling is developed and the new results 
presented. Model transition formulas relating the ARMA 
model to the MA and AR models are developed and the input 
signal requirements of the modeling process explored. It 
is shown that after suitable modification, the Levinson 
algorithm can be applied to solve the ARMA modeling problem 
recursively in order and lattice solution methods are also 
developed for both a batch processing and an adaptive model 
solution. The results of experimental simulations of both 
of these modeling solution methods are presented and dis- 
cussed, and comparisons are made with conventional means of 
ARMA modeling using the equation error formulation. Finally 
it is also shown that the lattice solution methods can be 
generalized to solve for the multichannel ARMA model with 
arbitrary numbers of inputs and outputs. 

A. LINEAR ARMA MODELING AND ITS RELATION TO AR AND MA 

MODELING 

The ARMA model for linear systems assumes the current 
value of the output of the system is given by a weighted 
combination of present and past values of the input and 
past values of the output. In terms of the discussions of 
Chapter I, Fao is assumed to be zero and Pao and Foo take 
on the following forms 

М 


FigLu (OO ] = 2, a(n) u(k-n) (3.1a) 


ПА 





M 


F,yLy(k-=1)1 = p b(n) y(k-n) (3.15) 


leading to a transfer function representation for the system 





given by 
M 
a(n)z ® 
a А E CN 
TS biz 
m= 


A number of methods exist for finding the model coef- 
ficients {a(n)} and {b(n)}. Ав stated in Chapter I, a MMSE 
solution via the direct form modeling approach requires the 
solution of a system of highly nonlinear equations and in 
general is untractable. An alternative is to first obtain 
an estimate of the denominator polynomial B(z) by some 
means such as AR modeling and then using this in the system 
shown in Figure 3,1, estimate the numerator polynomial A(z) 
by setting its coefficients to minimize the mean square 
value of the error. This method was first explored by 
Shanks. [Ref. 49] 

Another alternative is to apply the Prony method [Refs. 
8, 52 and 56] derived in Appendix E which obtains the model 
parameters by matching the impulse response of the system 
and model over the first N+M+l sample intervals. Both of 
these techniques share a common characteristic. They both 
start by independently estimating the denominator coefficients 


(or model poles) and then, given this estimate, solve for 
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Figure 3.1. Shanks method for ARMA modeling 
m@emmumerator coefficients (or model zeros). This is in- 


tuitively unappealing in that one would expect these two 
estimation problems to be more closely coupled with the 
zero estimates also affecting the estimates of the poles. 
The application of the equation error formulation to the 
ARMA modeling problem permits simultaneous estimation of the 
model zeros and poles and was first used by Kalman [Ref. 23] 
in work on self-optimizing control systems. The prediction 
error form of the model is considered here where Р. 15 set 


0 


me, zero while 


M 

Fj tutk)] = 2 an) ukon) (5 па) 
N 

FjLyOO = y(k)- 9  b(n) y(k-n) (3.2b) 


n=1 





The analysis model is depicted in Figure 3.2 where 
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Щи) = = a(n)z P? Ш) 
n=0 
and 
N 
ВН > sz" (з Б) 
п= 1 


апа these polynomials are the estimates of the system trans- 


mere rtunetion numerator and denominator, 
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Figure 3.2. The equation error formulation for 
ARMA modeling. 





The expression for the model error can be written as 


b 
T = 


| 
1 ‚ 43 ] A (З.ца) 


e (k) = y(k) ~ [y(k) 


E YCK) and b are defined as in equation (2.7) апа 


u (k) SEAS ө: a бе 


а? = [a(0) «++ acm]? (3.40) 


This results in an expression for the mean square error 
which is a quadratic function of the model coefficients, 


enn the MMSE solution for those coefficients given by 


| 1 
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ES (vu = -уу 
BE... ------- e 2 Je. (3.252) 
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Ec uu Ui ws 
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model Transition Relationships 
Comparing equations (3.4) and (3,5) with their AR and 
MA counterparts, it is clear that the ARMA model provides a 
generalization of these other models since they can be 
obtained from the ARMA model by assuming that either a or b 


is zero. Consequently it is susceptible to the same type oi 





bias introduced in the AR and MA models by the presence of 
additive noise on either the system input or output signals. 
To develop the relationships between these models further, 
consider the inversion of the correlation matrix in equation 
(3.5a) in terms of its component matricies. Since the left 
and right inverses of a nonsingular square matrix are the 


same, either 
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can be used to find the required inverse in partitioned form. 


solving for the right inverse of equation (3.6b) yields 


E -1 -1 
А 3 mn ED LA qc R T ) (3.7a) 
pU cu y 
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Solving for the left inverse of equations (3.6a) gives iden- 
NE "results for A and D while equivalent but different 


Eus are obtained for B and C given by 


E тї -1 -1 
VUL CD S U y im uu 
м 
аи -1 -1 
ES В. вв Дов. в (3.8b) 
Е Е = ye я Ey 
м yu uy 


Memms equations (3.7) in equation (3.5a), the solutions for 


the ARMA coefficient vectors are given by 


-1 -1 
b = (R, -R R R ) D 
ae + г 
Yu + иу уу 
cd -1 -1 
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yu а Ц yu у 
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Mae Macrix inversion lemma [Refs. 11 and 19] which states 
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Bor monssingular square matrıces E, G and E+FGH, can be 


used in equation (3.9) to rewrite them as 


-1 1 -1 -1 
D= R (R - R R R ) 
ER A V — A з, о ae 
e t. teen 
-1 
ER NE E ane 1 Сала) 
RERO. 
+ CN -1 zd Е 
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The all zero and all pole model solutions of corresponding 


orders however are 


г and b = В г 
= Y 


where subscripts are used to distinguish these solutions 
from their counterparts in the ARMA zero pole model. 


ЙӘШ equations (3.11) it follows that 


-1 -1 -1 
b = b MENU (R - В RR ) 
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Following a similar development, the left inverse relation- 


ships of equations(3.8) can be used to write 


-1 -1 
b = (R - Е Е R ) [> -R ce] co Же 
УР 0 “yu EE “uy —уу en ¿HA 
a. = (В ES ESTER ee 102712 
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Equations (3.12) are termed the "Zero Pole Model Transi- 
tion Formulas" and specify the relationships between the 
various models. It is interesting to note that equations 
(3.12a) and (3.12b) take the form of a linear observer with 
a new estimate of the solution given by the old estimate 
plus a gain times an error term. To gain some insight into 
the functioning of these formulas, consider the form of N 
and Ry ™ for the linear system described by the transfer 
Memetion of equation (3.1c) 
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N 
ЦИК в) 5 REOS Ry min) + ae) Ro frien) 
ізі L=0 


(8.13a) 
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Assuming that N=N and M=M, and writing equation (3.13a) for 


ВО --1 and equation (3.13b) for -M«n«0 results in 


quas wp aum 


г ЭЕ) 


These constraints on the system input and output auto and 
cross correlation coefficients are the ARMA modeling equa- 
tions of (3.5a) with the model coefficients replaced by the 
system parameters. In AR modeling, Dap is set to satisfy 

the constraints of equation (3.14a) with the assumption that 
a” is zero, The error term in the model transition formula 
(3.12a) then checks this solution to see if it also satisfies 


қ + 
E Econstraints of equation (3.14b) still assuming that a 


ШШ Zero. 


СЕОБЕ Е и Dap ea (15) 


81 





ТЕ this error is zero and the constraints of equation (3.14b) 


are satisfied, equation (3.12a) sets Бур = Dap and equation 


ай Е 
ПГ а to zero. ТЕ however the error is nonzero, 


(3.12a) adjusts Dap 


+ 
Бор ала ЕВЕ provides a Nonzero a7p* 


(3.12a) and (3.12d) are complementary, specifying the ARMA 


ШИЕ ШоОрогтісі to the error to obtain 


Thus equations 


or zero pole model solution of order M over N when given 
the N-th order AR or all pole model solution. 

In like manner, equations (3.12b) and (3.12c) give the 
ARMA model solution of order M over N when given the M-th 
order MA or all zero model solution. The all zero solution 
is obtained from equation (3.14b) assuming that 6 is zero. 
Equation (3.12b) checks this solution against the con- 
straints imposed by equation (3.l4a) with the same assump- 
tions and adjusts ae appropriately to determine 414, 
Equation (3.12c) then sets Бур» completing the zero pole 
model solution. 

Modeling Input Signal Requirements 

Another aspect of the modeling problem that must be 
considered is that of system identifiability. If, from the 
available measurements of signals, a model can be obtained 
that accurately represents the system's operation, the 
system is considered identifiable. The two issues that 
arise therefore are the measurement requirements (which 
signals must be measured) and requirements on the input 
signal used to excite the system during the modeling 


process. In the equation error formulation of the ARMA 
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Edel. both existing signals (system input and output) must 
be observed (or at least a knowledge of their auto and cross 
correlation functions must be available). Most discussions 
of input signal requirements for identifiability simply 
state that the system can be id&t ified if the input signal 
eeu ficiently rich, persistent or exciting. eg [Ref. 19] 
To explore the question of input signal requirements fur- 
ther, consider the mean square equation error cost function 
being minimized. Assuming that the equation error signal 
is ergodic and has finite energy its mean square value is 
obtained via time averaging as 


со 


E, = efe(k)*} = % em? Ga 


n=-0 


Applying Parsevals relation this becomes 


со T 
Е, = X eam = $ f Е(е??)Е# (е?) de (3.17) 
n = —0 ~ Т 
where 9=uT and * indicates the complex conjugate. The 


equation error is represented in the transform domain as 
ЕШ) - FBOZOHCZ) - ACZz)]IUCZ) @.19) 


and using this in equation (3.17), the cost function becomes 
TT 


| | лы” 4002 
ВЕ қ SG E | 85 (3.19) 
E 





showing that the power spectrum of the input signal acts as 
a frequency dependent weighting function on a transfer 
function error term. Therefore to identify the system 
equally well at all frequencies, the input must have a 

flat spectrum as will be the case for a white noise input 
or an impulse function input. Otherwise, the model trans- 
fer function will only be matched to the system transfer 
function over the range of frequencies where the input 
signal has significant power. 

As an example consider the equation error ARMA 
model for a fourth order system driven by a single sine 
Wave input at a frequency of m/3. According to equation 
meg), the model transfer function will only be required 
to match the system at this single frequency, and to accom- 
plish this, only a first order model is needed. Any in- 
crease in model order above first order therefore should 
have no effect, Figure 3.3a shows a comparison of the mag- 
ide spectrum of the fourth order system and its first 
order ARMA model obtained using the sinusoidal input and 
as anticipated they match at the frequency 1/3 (coinciden- 
tally they also match at one other frequency as well). 
Figure 3.3b shows the same comparison but with a fourth 
order ARMA model. It is clear that increasing the model 
order failed to improve its accuracy and that the model 
accurately represents the system only at the frequency of 


ШОС Signal and, by coincidence, at one other frequency. 
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3.3. ARMA models for a system driven by 
a Single sinusoid. 
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ШОШО be moted that this type of analysis to 
determine input signal requirements could be applied to the 


AR and MA models as well resulting in the same conclusions. 


Pee RECURSIVE IN ORDER SOLUTION FOR THE ARMA MODEL 

Since the equation error formulation for the ARMA model 
memaeceneralization of similar formulations for AR and MA 
modeling, it is reasonable to assume that a Levinson-type 
algorithm could be devised to obtain the ARMA solution 
recursively in order, and that from that algorithm, lattice 
filter methods applicable to the ARMA modeling problem could 
be derived. Attempts to develop such an algorithm directly 
for the ARMA modeling equation (3.5a), however, fail to 
provide useful results. The first problem that arises is 
in deciding which model order to make recursive; the order 
of the numerator polynomial, the order of the denominator 
polynomial or both. If it is assumed that the numerator 
and denominator are of equal orders (M=N), the ARMA modeling 
equations become as shown in equation (3.20). However, 
efforts to develop a Levinson-type algorithm for this system 
of equations, where the numerator and denominator polynomials 
of the model are incremented simultaneously to obtain re- 
cursive in order solutions, are still frustrated by the 
presence of the (N*1)-st row and column in equation (3.20). 
This arises because the numerator coefficient vector а” 15 
Ber) -vector while the denominator coefficient vector b 


EN vector. If it is further assumed that the coefficient 
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а(0) is known in advance or can be estimated in some other 
ion, equation (3.20) for the solution for the remaining 


2N coefficients becomes as written in equation (3.21). 


в (0) EUER EEN E (0) ВСЕ Е) R (1) R (1) 
УУ ҮҮ к yu уу үч 
Р 7 
| 
Вт)... к (0) |R_ (N-1) ,.. к (о) b (N) в (М) R (N) 
yy yu yu = -a(0)! yu 
par (0) к (1-М) |К (0) ТТІ а) RL) к (1) 
цу u | uu uu uy uu 
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y uy о uu uy | uu 
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The coefficient a(0) essentially has the role of a gain for 
the model, and a method for estimating it after all the other 
coefficients “are obtained (as was done in AR modeling) will 


be discussed later. Equation (3,21) can be written as 


R R b р р | 
=yy  =yu| |= ғ, yu 
p 22) 
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Now consider the form of a two channel autoregressive 


model where the two input channels are y(k) and u(k); that is 


y (k) 


x CK) 


x, (kK) =U CK) 


Using equations (A.4a) and (A.7b) the two channel AR modeling 


equations are 


іі 514 |4 4 уу -yu 
: C) 


Ruy Rau 921 220 Puy Suu 


ا 


E omparing eqiatlions (3.22) and (3.23) it is clear that 
with the exception of the gain term a(0), the ARMA modeling 
solution can be obtained from the two channel AR solution 

of (3.23). Furthermore equation (3.23) can be solved inde- 
pendently of the gain term via the Levinson algorithm as 
shown in Appendix A or via the multichannel lattice methods 
developed in the previous chapter. Then all that remains 

to complete the solution for the ARMA model is to estimate 
the gain term a(0) and solve for the other model coefficients 


using 


E00 Ето 
А a, 20 


-а(0) 


Jg 
A, 
A, 


a 22 
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Brom this it follows that the transfer functions A(z) and 
B(z) for the ARMA prediction error analysis model can be 
related to the twochannel AR prediction error matrix poly- 


nomial transfer function by 


Erz) 1 
= [I-Dez)] 25) 


-А(2) -а(0) 


It is also easy to relate the ARMA equation error to the 
two channel AR prediction error vector. Using equations 
(A.3) and (A.5), the two channel AR prediction error vector 


memwritten as 


T 3 
y(k) | а а 
e(k)? = БМТ |17707) (526) 
i" EO 
= 
= се, (к) е (21 
Defining 
il 
ф = (3,27a) 
К -а(0) 
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and postmultiplying equation (3.26) by Ed usus (3. 21)» xt 
ша Понс that 


( 


eco! y = y(k) - alO)uck) - CyGoli ucol] (3.275) 


11 
] 


} 


But from equation (3,4a) this is exactly the ARMA model 


equation error so that 

en(k) = e(k)” у 276) 
and 

с (е (КТ) = фр | (3.27а) 


EP rs the forward prediction error covariance matrix 

for the 2 channel AR model. Equation (3.27d) also provides 

a means of estimating the gain term a(0) after the two channel 
AR solution has been obtained by setting it to minimize the 


mean square value of equation error resulting in 


e(e, OO e, O0] 


5 (28) 
e le, Ck) } 


a0) = 


and completing the ARMA model solution. 
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The portion of the ARMA model solution in equation (3.24) 
ЕП Dy the two channel AR solution can be found recursively 
in order using either the Levinson algorithm or the lattice 
filter techniques. If the desired ARMA model order is not 
known in advance, a model gain term a(0) can be estimated 
for each order two channel AR solution to find the ARMA 
model of corresponding order along with its MSE. In this 
fashion, the entire family of ARMA models for the system from 
order zero to order N, along with their mean square errors 
are obtained, If, on the other hand, the desired ARMA model 
order is known apriori, the gain term need not be calculated 
at each stage. Only one gain term must be calculated to 
obtain the ARMA solution after the appropriate order two 
channel AR solution has been found. 

It has already been shown that to fully identify the 
system using the equation error ARMA formulation, the input 
Signal must have a flat spectrum as in the case of white | 
noise. When white noise is used as the system input u(k), 
simplifications emerge in the solution of the two channel 
AR model via the Levinson algorithm or lattice methods. 


р Е . 2 
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where h(n) is the sampled impulse response of the system. 


Bonsıider equation (A,i6c) for 


-1 


R8 (0) TT OOD в С) 0 
QD yy yu yy 
К я (3.29с) 
(0) FCR 0) КО. 0 
uy uu uy 
: CL) (1) Р 
fees shows that Кү» апа Ко? are zero and furthermore since 
n) and p(n) пешие Гор all n 1t 15 seen that 
—yu uu 
m) np 2 


for all n as well. This can readily be understood by con- 
Sidering the role of these two coefficients at each stage 


in the AR prediction of y(k) and u(k). k and k are the 


Ша 22 
coefficients used in trying to predict u(k) from past values 
of y(k) and u(k), and when u(k) is a white sequence it 
cannot be predicted, forcing these coefficients to be zero. 
No such simplifications occur in the backward prediction 
problem (and therefore in the K matrices) since even for 


a white u(k), a backward prediction of u(k-n) from subse- 


quent values of u and y is possible, This is because in 








general, a linear dependence of y(k) upon past and present 
values of u(k) can occur (and certainly will occur when the 
relationship between y(k) and u(k) is described by an ARMA 
model). 

As a result of these simplifications, it is seen from 
equation (A.19a) that the polynomials 4,202) апа а,,(2) аге 
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zero when u(k) is white and the ARMA model is given by 


Bez) 1-d,,(2) 0 I 
= 03230) 
-А(2) E | -а(0) 
In this special case it follows that 
B(z) = det[I -J) (z)] qoem 


and therefore stability of the ARMA model and of the two 
channel AR model are equivalent. In general,however, no 
Such connection exists for arbitrary input signals. Further- 
more, even when the system input u(k) is white, solutions 
Kor Ky and Ко will not in general be exactly zero since 
the required correlations are usually not known and must be 
estimated. 

This development showing that the ARMA model solution 
can be obtained from a two channel autoregressive model de- 
pends on two assumptions: 


1) The numerator and denominator polynomial orders in 


the model transfer function are assumed to be the 


ац 





same and are incremented simultaneously to build up 
the desired solution recursively in order; 

2) It is assumed that the coefficient of z to the zero 
power (a(0)) in the numerator polynomial of the model is 
either known in advance, or that another means of 
estimating it can be found so that it need not be 
estimated directly in the modeling equations of 
0209. 

The second assumption causes no concern since the two channel 
autoregressive solution is obtained independently of a(0), 
and given that solution, it has been shown that a(0) can 
indeed be estimated in another fashion in equation (3.28). 
The first assumption however, warrants further consid- 
eration since it seems somewhat restrictive (at first) to 
require that the numerator and denominator polynomials of 
the model have the same order when in fact, the system being 
modeled may have different order numerator and denominator 
polynomials. To see that this assumption is not restrictive 
in general, consider what is occuring as the model is built 
up recursively in order. At each model order n, the pro- 
cedure finds the best n-th order model (with n zeros andn 
poles) in a minimum mean square equation error sense. 
Making the model numerator and denominator orders different 
Mereequivalentiy, forcing some of the coefficients to zero 
in the model where the orders are the same) places a priori 
constraints on the model, forcing some of the poles or 


zeros to the origin in the z-plane, rather than allowing the 





model to place them at will to minimize the cost function. 
As an example, consider the process of obtaining an ARMA 


model for a system given by 


2 -n 
an) z 


H(z) = 1509. 
je b(n) 2Z" 
n=1 


where two of the systems four zeros actually occur at the 
origin of the z plane. Constraining any of the model zeros 
to the origin at orders one, two or three will result ina 
model with higher cost (MSE) than if they were not con- 
strained. Even at order three, a model without constraints 
Can be expected to use the "extra" zero to help in approxi- 
mating the effects of the system's fourth pole as shown in 
equation (2.13) yielding a lower cost and more accurate 
model than would result if one zero were forced to the origin. 
Only at order four are such constraints reasonable but even 
then, they are not necessary since the modeling procedure 
itself should recognize that the best fourth order model 

Will have two zeros at z=0, 

Therefore, it is seen that assuming equal orders for the 
model numerator and denominator is entirely reasonable as a 
general approach in obtaining MMSE models for unknown sys- 
tems or even reduced order models for known systems. When it 
1s known in advance that the best MMSE model for a system 


ШИ 2егов ас 2-0, ітпросіпс such a constraint on the model can 
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reduce the computational complexity of obtaining the solution, 
but even here the assumption of equal orders is not restric- 


уе. 


РФ ИА ТАТТІСЕ ҒОКМ АҚМА MODELING 

In chapter two it was shown that the lattice structure 
of Figure 2,7 could be applied to solve the multichannel 
AR modeling problem in terms of the reflection coefficient 
matrices given by equations (2.49). Since the ARMA model 
solution can be obtained from a two channel AR solution with 
y(k) and u(k) as the input channels, only the structure 
described by equation (3.27c) need be added to a two channel 
AR lattice to obtain the lattice form of the ARMA analysis 
model. It is interesting to look at the exact structure of 
this lattice model as shown in Figure 3.4 for a second order 
case. This structure is seen as a lattice interconnection 
of two single channel AR lattices operating on the input 
Signals y(k) and u(k). The coefficients on the main diag- 
Bnals of the K and K matrices specify the single channel 
lattices while the off diagonal elements specify the inter- 
connections, (This will also be the case for extensions 
to lattices with any number of channels.) 

The ARMA synthesis model implementing the transfer func- 
meom A(zZ)/B(z) can also be put in lattice form. The for- 
ward prediction in the two channel AR analysis lattice is 


described by 
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БЕЙЕТ оп (3.33) along with the equation of e IA апа 
the equations for the backward prediction (2.456) describe 


the structure shown in Figure 3.5 for a second order case. 


(2) 


To provide the required input at m , recognize from 


шащастоп (3.27с) that 


y 02) (2) 


е (К = е„(К) + а(0)е (к) (3.Зца) 
у 0 u 

If the ARMA model is an accurate representation of the 
system, the equation error eg CK) will be quite small (ideally 
zero) so that in general for the N-th order case 
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ern) 


This is indicated by the dashed feedback path in Figure 3.5. 
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ЕРМА MODEL SIMULATIONS; BATCH PROCESSING 

ARMA modeling procedures for linear systems using both 
the lattice filter method and a brute force matrix inver- 
sion with equation (3.20) have been implemented and the 
results of these two approaches have been compared in over 
a thousand model simulations of more than thirty different 
linear systems. The experimental results which follow are 
a representative sampling of these simulations. 

In the lattice filter method, equations (2.49) were used 
to calculate the forward and backward reflection coefficient 
matrices, with time averages over a specified interval used 


(0). 50) 


Memestimate the required correlations in P and 


Bn Eun «N-l. Equations (A.15), (A.16a) and (A.16b) 
were used to obtain the two channel AR model coefficients 
from the K and K matrices and with the gain calculated 

in equation (3.28), equation (3.24) was used to obtain the 
desired ARMA model coefficients. Equations (A.18) were 

used to update the forward and backward prediction error 
covariance matrices from one lattice stage to the next. 
Figure 3.6a provides a flow diagram of the procedure. 

In the brute force matrix inversion method with equation 
(3.20) time averaging was again employed to estimate the 
required correlation coefficients. A rectangular window 
Was applied to the data, however, to retain the even symmetry 
of the autocorrelation function in these estimates. The 


ARMA model coefficients were then obtained using a general 


purpose library subroutine (which employed gaussian 


LOL 





elimination) to solve equation (3.20). Figure 3.6b provides 
a flow diagram of this procedure. 

In both cases zero mean, unit variance gaussian white 
noise was used as the system input. In the simulation re- 
sults that follow, a description of each system discussed 
eeen rer function coefficients, zero locations and pole 
locations) is listed in tabular form and root locations in 
the z plane as well as transfer function magnitudes are 
plotted for the system, and various models obtained for it. 
In each case, models were obtained for averaging intervals 
2200, 500, 1000, 2000 and 5000 data points. Only the 
results for the two extremes of 200 and 4000 points are 
included here. 

The first system considered hasasecond order numerator 
in Z inverse and a fourth order denominator, and its char- 
acteristics are listed in Table 3.1. Figure 3.7 shows a 
comparison of the root locations and transfer function 
magnitude of this system with those of fourth order lattice 
filter and brute force models obtained when correlations 
were estimated by averaging over 200 samples. Figure 3.8 
provides the same comparison for a longer averaging interval 
of 4000 samples of data. While both methods perform com- 
parably with the longer averaging interval, the lattice 
Ed produces a far more accurate model with the short 
Berosıng interval. It is also interesting to compare the 
performance of the two methods when the system is overmodeled; 


that is, when the model order is increased beyond that of 
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Figure 3.6a. Flow diagram of batch processing 
lattice solution for all ARMA 
models orders 1 to N. 
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ABLE 3.1 


SYSTEM A 
RAN ИБОАСТІСМ COSFFICIENTS 
NUMERATOR D=NOMINETOR 
А(С)= 0.25000 
А (1)= 0.35000 8(1)* 1.14000 
Al2)= 0.24500 8(2) = -1.45490 
A(2)= 0.0 B(3)= 0.88490 
А (4) = 0.0 8(4)= -0.40745 
ROOT LOCATIONS 
ZEROS POLES 
RE IM RE IM 
-0.7СС00 С. ТОООС 0.50020 0.50000 
-0.7СС00 -0.70000 0.500096 0.50006 
С.С 0.0 0.07090 0. 90000 
С.0 0.0 0.07009 -0.90000 
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ЕЕ узсет. Figures 3.9 and 3.10 show such a comparison 
when sixth order models are obtained for this fourth order 
system. Ideally, the extra zeros and poles should lie at 
the origin in the z plane making the additional transfer 
function coefficients zero. Figure 3.9 shows that for a 

200 point averaging interval, the lattice locates the extra 
roots in the vicinity of the origin while the brute force 
model does not. When the averaging interval is increased 

to 4000 points, the extra roots of the lattice model move 

in toward the origin while those of the brute force model 

do not. Instead, a zero pole cancellation at some arbitrary 
location occurs in the brute force model. In all cases 
investigated during this effort, the lattice method clus- 
tered the extra roots in the vicinity of the origin and as 
the averaging interval was increased to take in more data, 
these roots were consistently moved in closer to the origin. 
This property is further evidenced by the plot of mean square 
equation error as a function of model order shown in Figure 
ШИ Гог a 200 point averaging interval. The MSE for the 
lattice model flattens out at order four indicating that 
further increases in model order fail to increase its 
accuracy. Meanwhile the MSE for the brute force model con- 
tinues to decrease beyond fourth order as it uses the addi- 
tional roots to reduce modeling errors caused by innaccuracies 


Ши Ее fourth order model. 
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Figure 3,11. Mean square value of equation error 
(as a percentage of the mean square 
value of system output) vs model order 
for lattice and brute force models of 
system A with 200 point averages. 
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To investigate the ability of these modeling methods to 
distinguish roots located near one another in the z plane, 

a pair of zeros were added to the previous system in close 
Brozimity to one of the pole pairs. The characteristics of 
NEEESvstem are listed in Table 3.2. Figures 3.12 and 3.13 
show the lattice method and brute force modeling results for 
200 and 4000 point averaging intervals. With data over only 
200 sampling instants, neither method is able to accurately 
model the effects of the adjacent roots. When the averaging 
interval is increased, the lattice correctly models the 
Plant while the brute force method does not, and even results 
in an unstable model. (Figure 3.13b comparing the transfer 
function magnitudes has been plotted in spite of the model 
instability.) 

To investigate the ability of these zero pole modeling 
methods to model systems that are actually all zero or all 
pole, the systems listed in Tables 3.3 and 3.4 were used. 
The modeling results are shown in Figures 3.14, 3.15, 3.16 
mma 3.17. 

The conclusions drawn from the results of this experi- 
mental study are as follows: 

ie tor short data lengths, the lattice filter method 

provides more accurate models than does the brute 
force modeling method. 

2) When the system is overmodeled using the lattice 

method, the excess roots are.clustered about the 
origin in the z plane and as the averaging interval 


ls increased and more data taken in, these roots move 
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IBIDEM 3.2 


So TEM B 
INS ERA UINET LEN COEFFICIENTS 
NUMERATOR DENOMINATOR 
А(Сіс 1.С0000 
А(1)= 1.34000 B(1)= 1.14000 
А( 2): 1.79940 B(2)= -1.45490 
A(2)= 1.20600 2(3)= (0.88490 
A(4)= 0.88533 8(4іс -0.49745 
ROCT LOCATIONS 
ZEROS PCLES 
RE IM RE UM 
-C.7CCCC 0.7000C 0.50000 0.50006 
ЕС. ІСС0ОС -0. 700CC 0.50000 -0.50000 
0.02000 0.5500С 0.07000 С.5000С 
C.C3CCO 2:09 9500€ 0.07000 -0.90000 
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ПАВЕЕ 3.3 


SYSTEM C 
TRANSFER TFUNCTICN COEFFICIENTS 
NUMERATOR DENCMINZTOR 
А(С)= 1.60600 
A(1)= -1.50000 8(1%- С.С 
Al2)= C.62500 B(2)= 0.С 
ROOT LOCATIONS 
ZEROS FCLES 
RE IM ВЕ IM 
0.75000 0.25000 0.0 0.0 
С. ТЕССС -0. 2500C 0.0 0.0 
TABLE 3.4 
SYSTEM D 
TRANSFER FUNCTICN COEFFICIENTS 
NUMERATOR DENCHINETOR 
А(С)= 1.560600 
A(1)= 0.0 B(1)= 1.50000 
A(2)= С.О 8(2)- -0.€2500 
ROOT LOCATIONS 
ZEROS POLES 
RE IM КЕ IM 
0.0 0.0 0.75000 0.25000 
0.0 0.0 0.75090 -0.2500C0 
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consistently in toward the origin. This also forces 
the excess transfer function coefficients and re- 
flection coefficients to be very small (ideally zero), 
clearly indicating that the model order is higher 
than necessary. The brute force method,on the other 
hand, scatters the excess roots throughout the z 

plane and produces cancellations of the excess Zeros 
and poles. This results in nonzero values for the 
excess transfer function coefficients. 

3) The MSE as a function of model order is generally 
better behaved for the lattice method than for the 
brute force method, decreasing rapidly until the 
correct model order is reached and then failing to 
decrease substantially as the model order is increased 
Further. 

Two qualitative explanations for the improved performance 

of the lattice filter modeling method are offered. The 

first is that the data 1s not windowed in the lattice method 
while a window that is nonzero only over the span of the 
averaging interval must be used in the brute force method to 
maintain even symmetry in the autocorrelation function esti- 
mates. The effects of this difference between the two modeling 
methods should be most noticeable for short data lengths, and 
become less evident as the length of the averaging interval 
is increased. The second possible cause for the lattice 
method's better performance is that the actual output se- 
quences of the n-th order lattice are used to calculate the 


coefficients at the (n*1)-st order stage. In this manner, 








modeling errors that have occurred in the n-th order lattice 
can be compensated for to some extent in the (n*1)-st order 
stage. No similar phenomenon is evident in the brute force 
modeling approach. Consider the differences between the 
Levinson algorithm (which is equivalent to the matrix inver- 
sion method) of equations (A.16) and the lattice method given 
by equations (2.49). Both methods calculate the corrections 
that must be made to the n-th order model to obtain the 
(ntl)-st order model in terms of кош апа wu) and 


theoretically, 
(n), zn) D T (п) pfo 
Бс k-D |-R (nt) -p ) 


When correlations are estimated by averaging over finite 
intervals however this equality will not in general be 
Satisfied making the two methods different. The Levinson 
algorithm will estimate the correction terms to be added to 
the optimum n-th order model while the lattice method esti- 
mates the correction terms to be added to the estimated n-th 
order model actually obtained. 

The improved performance of the lattice method is not 
achieved without cost, however. The method is made compu- 
tationally expensive by the need to store the system input 
and output sequences and the lattice prediction error se- 
quences and pass them through successive stages of the 
lattice as it is built up in order. The computational com- 


plexity of the two methods is compared in Table 3.5. 
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Mable 3,5. 


Number of 
Correlation 
Estimates 
Required 


Matrix 
Inversions 


Data Storage 
Requirements 


Computations 
M Pass Data 
rough The 
Filter 





Computational requirements for batch 
processing ARMA modeling of an N-th 
order system using P samples of the 
system input and output. 


uUN+3 averaged over 
P data points 


2N-dimension 2 
l-dimension 1 


1 - dimension 2N*1 


2P samples 


8NP multiplica- 
tions 
UNP additions 
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ADAPTIVE LATTICE ARMA MODELING 

In addition to the batch processing method described in 
the previous section, the lattice ARMA analysis model can be 
implemented adaptively as well. The adaptive lattice solu- 
tion for the multichannel AR lattice, which solves most of 
the ARMA modeling problem, has already been described in 
chapter II. To make the lattice ARMA model adaptive, only 
an adaptive solution for the gain term a(0) need be added. 
To avoid ambiguity, the time varying adaptive estimate of 
this term at time k is denoted by ay fk). Applying an LMS 
adaptive algorithm it follows that 


ag(k*1) - ag Ck) = на V(K) COS) 


0 


and using equation (3.27d) to form an instantaneous estimate 


of the gradient yields 


КОШОК a2 e (k)le Ck)-a.CKkJe. Ck) ] (3.36) 
u y 0 u 
= -2 e, (Kk) eg CK) 
so that 
ay (k+l) = ag (X) Б 205 е (К) eg OX) З) 
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Here it is clear that 
1) е СК) is analagous to the error signal. 
2) е (к) is analagous to the input signal. 
B э is analagous to the desired signal. 


ЕР Спас for stability Hg must satisfy 


Е Уе (3233) 


Once again, however, the mean square value of e, 09 may vary 
from stage to stage and over time as the two channel AR 
lattice adapts making it appropriate to apply relations 


lar to equations (2.58) and (2.59) 


Но (k) - (КУ (32344) 


ео (шту чы е (к) (3.39b) 


0 0 
where a is the normalized adaptive step size and the depen- 
dence on order is implicit (a superscript (n) could be used 


on a anda line ervor terms 1n equations (3.35) 


0 90 
through (3.39) to explicitly denote their dependence on the 
Brder of the solution), 

This adaptive lattice ARMA modeling scheme was implemen- 
ted and the results of its use in modeling system A described 


in Table 3.1 are presented here. A normalized adaptive step 


Size of a = .05 was used in the following simulations and 
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the results represent an ensemble average of one hundred 


trials. Unity variance white noise was used as the system 
input. A flow diagram of the procedure is shown in Figure 
FTB. 


Figure 3.19 shows a plot of the mean square equation 
error as a function of time for a fourth order model while 
1t adapts and Figure 3.20 shows the behavior of one term in 
(п) > —(n) 
11 ? and one term in the K matrix, Куу» 


Irse five lattice stages, 1 <n <5. These graphs of 


the K matrix, k at 


Ee Koi terms clearly show the successive stage by stage 
manner in which the lattice model adapts. Figure 3.21 shows 
a comparison of the transfer function magnitude and root 
locations of the system with those of the fourth order 
adaptive model after 1000 and 2000 iterations. Figure 3.22 
shows the same comparison in the overmodeled case when a 
meth Order model is obtained for this fourth order system. 
While these results show that the adaptive solution 
performs well and is a viable alternative to the batch pro- 
cessing solution, Figures 3.22e and 3.22d show that one 
advantageous characteristic of the batch solution has not 
carried over. The excess roots in the overmodeled case are 
ШЕР ГІРПТІу сіпстегес іп the vicinity of the origin indi- 
cating that the excess transfer function coefficients are not 
near zero. In examining Figure 3.20 it is also evident that 
the convergence of the reflection coefficients at the 
(5) —(5) 


( 
11 апа К+ 


Value of zero is quite slow. This relatively slow tracking 


overmodeled lattice stage (k ) towards their true 
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of overmodeled, zero valued parameters has been found to be 
a general characteristic of the adaptive lattice algorithm 
and has also been noted briefly by Morf. [Ref. 38] 

To understand the reason for this behavior, consider 
what occurs as the overmodeled fifth stage in the previous 
example adapts. Initially, before the coefficients of the 
first four lattice stages converge to their optimum values, 
the prediction error sequences out of the fourth stage are 
large and suboptimum. These signals provide incorrect in- 
puts to the fifth stage driving its parameter estimates to 
some values other than their optimum zero values. As the 
first four stages converge, the prediction error sequences 
ing into the fifth stage get small and since they drive 
the gradient estimates, convergence back toward zero is 
quite slow in these parameter estimates. 


The cost functions being minimized at the overmodeled 


(5) 5(5) 


И 5сасе are the trace of P CM given by 
T T Ш 
ро E E _ бз cC» 2 452 д Б «о» Bel 
(3.40a) 
T T T 
Бо E ES nd к nd oe ۴ jm ABL 


(3.40b) 


Applying the results of Appendix F, the parabolic surfaces 


defined by these cost functions are described by the eigen- 


(4) BD 


Eenes and eigenvectors of P and , the prediction error 


covariance matrices at the fourth stage. Consider the forward 
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Beedictions. The actual system output is given by 


ц ц 
= > 5) Ук) + 2, а(1) ч(к-1) Cul 
L=1 1=0 
and for a white input signal, the minimum errors in the 
fourth order two channel autoregressive predictions of y and 


a ane 


ey ІІ CP (3.42a) 
е") Ge SECS (3.426) 


This results in an optimal prediction error covariance matrix 


given by 


00) aco 


ВО = ШӘ COD (3.43) 
ще uu 


a(0) aL 


with eigenvalues of 0 and Шар)“. Also in the case of the 
system described by Table 3.1 where a(3) = a(4¥) = 0, it is 
seen from equation (3.41) that a perfect backward two channel 
AR prediction of y(k-4) is possible resulting in an optimal 


backward prediction error covariance matrix of 


(4) ` 


|| 
и 


(3.44) 
0 se 


тр е 








which also has a zero eigenvalue. 
Since the ellipses obtained by passing a plane through 
the parabolic cost surface have axes whose half lengths given 


(4) (4) 


Бу 1/УХ, ald Since Р and P each have one eigenvalue of 


52) 


Ш ЕО, ІС 15 seen that the parabolic bowls along which K 


ana KË? 


adapt, degenerate toward infinitely long parabolic 
troughs as the first four stages converge toward their opti- 
mum values. This is responsible for the slow convergence of 
the overmodeled parameters back toward zero. To avoid this 
problem, some means of detecting this degeneration of the 
cost surface and then reseting the appropriate parameters 


to zero must be found and this certainly provides an inter- 


esting area for future study. 


Е. A LATTICE APPROACH FOR MULTICHANNEL ARMA MODELING 

The lattice filter solution methods for the ARMA model 
can readily be generalized to multiple input multiple output 
ARMA models for linear systems. The equations for the 
multichannel ARMA model of the system shown in Figure 3.21 
are developed in Appendix G with the solution for the model 


coefficients given by equation (G.7) repeated here for 


Ponvenience. 
R R B К 
ES УТ" -- —Y Y 
= 9757 
R R А” R 
"uy ШІН г "Uy 
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Figure 3.21. A general multichannel ARMA system 


This 1s clearly a generalization of the single channel ARMA 
equation error solution given by equation (3.20) and just 
as in the single channel case some assumptions are required 
to apply a Levinson type algorithm, 

If it is assumed that all the Де are known or can 
be estimated in another manner, they can be incorporated 


into a matrix given by 





a,, 69) ay (0) 
Ҹо 
Ag = С ще) 
а GOD avec. а (0) 
91 9190 
al 
and equation (3.45) becomes 
М 50 |2 E Ryy Ryu 1 
Cai) 
Ем Eu | 4 ST TÉ Ге 











= Similar in form to the equations for the solution of 
а ©; + Qo) channel autoregression with input channels уук), 
EN (5092, ш (Оз eS им CK) so that the multichannel 
90 = Ro 
ARMA solution is related to the multichannel AR solution by 


- D (3.48) 


A eo 


The multichannel AR prediction error vector is given by 


Jt a 
YK) e (k) 
| 
е(к) 1 |----- - ра | 01 D = ----- (3.49) 
u(k) | Е e (X) 
and defining 
I 
р = (3.50) 
-А 


TS follows that 


e(k) ty Coz ca) 


Ж m | 
ZOD = uD Ag = [Y U] 


е С)“ 
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establishing the generalization of equation 3.27c relating 
Beenultichannel AR prediction error vector to the multi- 
channel ARMA error vector е (к), The ARMA prediction error 


covariance then is given by 
үр 0552) 


and the coefficient matrix An CAR DESE SEC Lo mınımize the 


Ei P- resulting in a solution given by 


0 


-1 


EE e (k) e G0! 
= =u =u 


0 ele, (k) e, 00) в.) 


completing the multichannel generalization of the single 
channel results. The portion of the solution given by the 
multichannel autoregression can be solved as before using 

the lattice methods in either batch or adaptive fashion. 

Then the matrix of gains can be obtained from equation (3.53) 
by batch processing or equations (3.37) and (3.39) can be 


generalized to yield an adaptive solution given by 





А„(к+1) = Ag(k) + 2 е (к) e coe а) 
- - Zu —0 
с СК) 
where 
2 _ 2 T 
с С) = (1-а) 24 0-1) + ае (К) е (К) (3,54b) 
-ц Zu 
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It is clear that the equations and methods deyeloped earlier 
for the single channel ARMA model are a special case of 


these results with Q, = Q, = I. 


TOE 


гуша >» = ср 
г 





IV. NONLINEAR SYSTEM MODELING 


The modeling of nonlinear systems is a far more complex 
problem than linear systems modeling. No attempt is made 
here to provide a comprehensive treatment of the problem. 
Rather, two specific models for systems comprised of the 
interconnection of linear and memoryless nonlinear subsys- 
tems are considered. Both of these models, the Volterra 
model and the new nonlinear ARMA model, are shown to be 
generalizations of the MA and ARMA modeling problems explored 
previously so that with appropriate modifications, the 
Levinson algorithm and lattice methods can be used to solve 


for the model parameters. 


A. VOLTERRA NONLINEAR MODELING 

The Volterra series model characterizes nonlinear sys- 
tems using a generalization of convolution where the system 
output is approximated as a summation (possibly infinite) 


Of Outputs of degree m systems. 


^ М 
Bep xo) (4. 1a) 


ме | 


This is shown pictorially in Figure 4.1. 


19,2 











Del) 
M-th degree 
system 


u(k) 





l-st degree 
SyS Cem 


Figure 4.1. The Volterra model for nonlinear systems 


A number of representations for these m-th degree systems 
are possible. The most commonly used representation is in 
terms of symmetric Volterra kernels and is given Бу 

со 


СЕ) = 2 s У a,(n]...n Ju(k-n,)...u(k-n_) (4.1) 
n= 9 n4-0 


and а. (п Па) is the m-th degree symmetric Volterra kernel. 


1° 
(Any permutation of the indicies results in the same value 
for this kernel giving a high degree of symmetry.) 

This model arises quite naturally for a linear system in 
cascade with a power series nonlinearity as shown in Figure 


4.2 for a quadratic nonlinearity where the output can be 


Dritten as 
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ШК) - » >» a(n, a(n, )u(k-n, )u(k-n,) 2 20 
n,=0 n520 
and 
a ,(n4n5) = a(n, )a(n,) (4.2b) 










со 


Ak n anul k-n) 
n=0 


ЕСК) (469 


Figure 4.2. А quadratic nonlinear system. 


The Volterra series model in this form has been widely dis- 
Ж ЕС in the literature Leg. Ref. 1, 6, 13, 14, 25, 26, 48 
and 54] since Wiener l|[Ref. 62] first applied it to systems 
analysis and modeling. ТЕ has the added benefit of treating 
linear systems as a special case of the model since the 
first degree kernel is exactly the convolutional represen- 
tation of the MA model. 

The number of kernels (M) required for an accurate model 
depends on the nature of the nonlinearity in the system and as 


long as the nonlinearity is soft, a relatively low degree 
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model will suffice, keeping the problem manageable in that 
regard. The primary difficulty associated with Volterra 
nonlinear modeling arises from the fact that it uses a 
nonrecursive MA representation for the linear portion of 
the system, requiring in general, an infinite memory as 
indicated by the upper limits on the summations in equation 
(4.1b). In practice, these summations need to be truncated 


as shown in equation (4.3) 


а ок 7... (п ) A43) 
5 m Ш m 


N 
ШЕ) - Dur 1 
n4 -0 n 


3 “= 


0 


but a large number of terms may still be required to 
accurately model the system. One method of solving for the 
model is to set the parameters (terms of the Volterra kernels) 
to obtain a minimum mean square equation error where the 
equation error is defined as the difference between the sys- 
tem output and y (k). 

To simplify the solution and reduce the number of para- 
meters that must be obtained, the symmetry of the kernels can 
be exploited by rewriting equation (4,1b) as 


y (xk) = » B uus 25 ере 0-07 
a т 


0 Ny =n, 


(4.4) 


where а, (п Dy) is the m-th degree triangular Volterra 


1: 
kernel. For a finite upper limit of N on the summations 


the m-th degree symmetric kernel will contain (N*1)" terms 


165,5 





| 
Bas only OT of them are distinct with the remainder de- 


termined by symmetry considerations. 

Still another representation for y, CK) uses the regular 
form of the Volterra kernel (this terminology has recently 
been Pe coduced ЕЕ, апсу апа РарРЬ [Веї<. 7 апа 35]). 


It 1s given by 


у (К) = 2 T 2; ME Зы: (кам =i) 
В. = h_=0 
1 m 
u(k-hi-...-hj) (4.5) 
where a (hy ..-h,) is the m-th degree regular Volterra kernel. 


With infinite upper limits on the summations, the symmetric, 
triangular and regular forms of the Volterra kernels are 
equivalent, however, when finite upper limits are used they 
cover the field of the model kernels in different ways. 
Because of its symmetry, it is reasonable to have equal 
upper limits on all the summations on the symmetric kernel 
as was done in equation (4.3). This is shown in Figure 4.3a 
for a second degree kernel. The equivalent triangular kernel 


is given by 


N 
ШЕСІС) - >= 


N 
m = a,(n]...n Ju(k=n,)...ulk=n, ) (4,68) 
n = 


en ned 


and it covers the region shown in Figure 4.3b for a second 


degree case. The corresponding regular form expansion, however, 
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requires variable upper ni omethe summations to cover 
the equivalent kernel space as shown in equation (4.6b), and 


Figure 4.3c for a second degree kernel. 


N — N-h, N-h, 
y_{k) = pe ec Y ОООО Du(k-h,)u(k-h, -h5) 
l 2 An 
(К-у...) CU. 6B) 


Thus it 1s seen that only half the field needs to be covered 
by the triangular and regular expansions to identify the 
kernel associated with the square field of the regular ex- 
pansion. When the regular form expansion is used with con- 
stant finite upper limits there is no inherent reason to 


make all the upper limits equal since there is no symmetry 


in the kernel. Considering the regular expansion 
m m 
у (К) = NECS ие) 
р. =0 h =0 
| m 
(4.7a) 
a triangular expansion given by 
PANIN 1o 
а "an 
y5Ck) z 2 E a, (ny ---n Jutk-=n,),..ufk=n,) 
ny = са 1 
СШ? 


Bemmequired to cover the corresponding field, For a quadratic 
expansion this is shown in Figures 4.4c and 4.4b, The equi- 


valent region in the symmetric field 1s shown in Figure 4.2a. 
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Therefore, identifying a rectangular kernel in the regular 
kernel space is equivalent to obtaining a symmetric kernel 
in the arrow shaped region of Figure 4.4a. 

Which type of expansion is more appropriate for a given 
system depends on the shape of the kernels for that system. 
For example, if a quadratic nonlinear system has a kernel 
with a relatively square shape in the symmetric kernel 
Space, a regular form expansion with constant upper summa- 
tion limits will have to estimate many zero valued terms 
and is inefficient. On the other hand if the system has a 
kernel similar to the arrow shaped region of Figure 4.4a 
in the symmetric space, a regular expansion will be efficient 
while a symmetric expansion over a larger field would be 
required with many zero valued parameters. Not enough is 
known about how to relate types of nonlinear systems with 
various shaped kernels however, and the question of kernel 
shape and the best type of expansion is not pursued further 
here. 

1. Lattice Filter Methods For Volterra Modeling 

Lattice filter methods derived earlier can be applied 
to obtain a minimum mean square equation error solution for 
the Volterra model if the regular form of the expansion is 
used. The Volterra model can be put in the form of a mul- 
tiple input single output MA model by defining a new family 
of signals as nonlinear combinations of delayed values of 


the system input u(k). 
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(k) =uCk) u(k-h,) u(k-n с ка 2222514) 


2 2 


h, ШЕ 


(4.8) 


For finite summations the regular form of the expansıon 


becomes 
ev №2 s 
y, = Dy КУ D» > a eae len 
Dc h,=0 h,=0 А m 


(+.9) 


ШО Сап Бе regarded as the sum of the outputs of a large 


number of linear filters whose inputs are given by u co 


h Be 
2 m 
Equation (4.9) is exactly the form of a 0. input, single 


output MA model where 


Q, = (N,,*1)(N. 


ЛЛ М. FI) ВО) 
B mm 


3 


Furthermore, if the same upper limit on the summations over 
hy 1s used in each of the various m-th degree systems, 
(МТ =М. 1 =. --=М 1), the overall Volterra model given by 
equations (4.la) and (4.9) is in a form suitable for solu- 
tion via the multichannel Levinson algorithm or the multi- 
channel MA lattice methods. The requirement to use the same 
upper limit on all summations over hy arises because the 
Levinson algorithm and lattice methods assume the same amount 
of memory in each channel of the model, 

As a specific example consider a second degree ex- 


Pansion where N =N, and N,,=N 


in, ege sth DERE 


ToL 





N N 


= 


1 
y(k) = X a,(hdu(k-h,) + >, КОЮ 
h, =0 0 2 

T 2 1 
(4.11а) 
ц. (к) = u(k)uk-h,) (4.11b) 
2 


Defining data and coefficient vectors for each channel given 


by 


+ _ Т 
цу (К) = Гы, (k)...u, (k-N,2] (4.12a) 
2 2 2 
г 3 Te 5) (N, hu) 1 (4.12b) 
ің, bee мр a 


and embedding them into single data and coefficient vectors 
written as 


T 


X o% s qe 


и 
га 
{= 
E 
- 
ІС 


E! | ШЕ 20) 


а, 
11 
=a 
f» 
| tv 


the equation for the Volterra model output becomes 


у(х) = Х* (к) 


¡a 


(4.12e) 


16 2 





which is clearly of the same form as equation (A.24a) and 


represents a N,*2 input, single output MA model. All the 


2 
nonlinearities in the Volterra model are external to this 

MA model, in the formation of its various input signals from 
the system input u(k). This model is illustrated in Figure 
E. 

It is interesting to consider what the recursive in 
order nature of the Levinson algorithm or lattice methods 
mean to the nonlinear Volterra problem. In building up the 
MA model solution recursively in order, the upper limit of 
the summations over hy 1s increased until the desired value 
1s reached. In terms of the Volterra model kernels, this 
means allowing each of the regular form kernels to grow in 
size in the hy dimension while holding their boundaries fixed 
at prespecified values in all the other dimensions. In the 
linear MA model, the kernel has only one dimension and there- 
fore the recursive in order solution eliminates any require- 
ment to prespecify its upper limit (the order of the model). 
For the higher degree nonlinear kernels, these: methods reduce 
by one the number of kernel boundaries that must be pre- 
specified for each kernel. Allowing the kernel to grow in 
any of the other M-1 dimensions (n, through h corresponds 


to adding additional channels to an existing lattice and 


simple methods of accomplishing this are not available. 


11515 





y5Ck) 


200) 





Figure 4,5. A second degree nonlinear Volterra 
model in multichannel MA form. 
ИШСЕ Л Кш гесе в ane Cocifi- 
clients" of the single input single 
output MA models shown. 
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В. NONLINEAR ARMA MODELING 

As was previously mentioned, the primary difficulty 
associated with the Volterra series model arises from the 
fact that it is a nonlinear generalization of the MA model 
and as such, a large number of terms may be required to 
accurately represent even a mildly nonlinear system. In 
linear system modeling this difficulty was remedied by 
using the more general ARMA model. It is reasonable to 
assume therefore that a nonlinear generalization of the 
ARMA model could remedy the problem in the nonlinear 
modeling case (for at least certain types of nonlinear 
systems). Such a generalization called the nonlinear ARMA 
model has recently been proposed. [Ref. 64] This model 
forms an estimate of the current value of the system output 


as follows: 


n.=0 


ШІ S^ s (nOu(k-n,) * 22 2 a (nn )ulk-n, )ulk-n,) 
S ШІ 8 u F = 1.2 a 2 
1 п-50 n,=0 


5 2, E ursi een) 
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+ DE b (m,)y(k-m,)* » 2, b (m,m Jy (k-m, )y(k-m,) 
_ cM JE = E Sl? i 2 


+ + cA = 
us уйк UP b, n. т ЗУСК m4)... y(k E 
т. =1 
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Ir + do clas 


n, =9 т. =1 n,=0 pos ml 


A DE A Dl kna ul k-n )y(lk-m,) 
m71 I pP l al Ë р ЩЕ 


КМ 4.13 
y( т ( ) 


The first three terms of equation (4.13) are a discrete 
Volterra expansion of the input signal u(k) and represent 

FP, luk) J im terms of the discussion of Chapter I. The 
second three terms are a discrete E EE expansion on the 
system output delayed one sample interval and represent 
Fogty(k-1)]. The final two terms are bivariate expansions 
of the system input and delayed output and represent 
Pluk), pao ENIS iS Té first and only model con- 
sidered where Pyle J 1S not assumed to be zero.) Equation 
(4.13) is clearly a nonlinear extension of the linear ARMA 
model contained in the first and fourth terms. As was the 
case in the Volterra model, the number of multiple summations 
required in (54.13) is dependent upon the nature of the non- 
linearity in the system being modeled. The upper limits on 
the summations in (4.13) have been written as infinity in- 
dicating a requirement for infinite memory or model order. 
As is discussed subsequently, however, the required model 
order (memory) may in fact be finite due to the nature of 
the system being modeled, thereby alleviating the difficulty 
encountered in the Volterra nonlinear model previously presented. 


EE 





The kernels of the input expansion and output expansion 
а) апа DD are symmetric since any permutation of the 
indicies results in the same value for the kernel. They 
have therefore been labeled with a subscript "s" and can 


also be written in triangular and regular form. 


D.» a (nj...n )u(k-n42...u(k-n ) 
o F р р 
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hy 50 ПРЕД НЕ al q | ^ 
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вит спе Ее triangular and regular forms in equation (4,15) 
the lower index on the summations has been shifted to zero. 

In the case of the bivariate expansion terms in equation 
(4.13) the kernels do not possess symmetry so that trian- 


gular expansions are not possible but regular form expansions 


are possible. 


2 пт о, 
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1 ptq 
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= p E. 21 со (ћу py dy doy). oy (lo... A) 
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Ш ptq 
u(k-h, 7... cho 42. u(k-h-. В ра) 
(4.16) 


Thus two regular forms are possible for the bivariate expan- 
Sions, Figures 4.6 and 4.7 illustrate the manner in which 


the regular form of the bivariate expansion covers the field 
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of model kernels for a second degree case for finite upper 
mits of Ny and М, оп the summations. It is interesting 
to note that because of a lack of symmetry in the original 
form of the bivariate expansion, the causal region in the 
regular form extends outside the quarter plane. 

As was the case with the Volterra model, it will be 
shown that in regular form, a minimum mean square equation 
error solution can be obtained for the nonlinear ARMA model 
coefficients using either the Levinson algorithm or the 
lattice methods. This will provide the nonlinear generali- 
zation of the results presented in Chapter III on linear 
ARMA modeling. Before developing this method of solution, 
however, the applicability of the nonlinear ARMA model to 
various types of systems and its memory requirements will 
be considered. 

1.  Identifiability Conditions and Memory Requirements 

In the previous chapter on linear ARMA modeling it 
was stated that there were two facets to the question of 
model identifiability; input signal requirements and mea- 
surement requirements, The question of measurement re- 
quirements was not discussed, however, since it was assumed 
that all signals (input and output) were observed. In the 
study of systems comprised of interconnected linear and 
nonlinear subsystems, however, various internal signals exist 
and the effect of either observing or not observing them on 
the modeling process must be explored. To do so, it will 


be assumed that the system under study can be put into the 


form of Figure 4.8 fulfilling the following equations. 


n 





х (К) = Ii уу + C,u(k) (4.173) 

xu CK) = Г. ур E C,u(k) (7D) 

O =F [x (x)] (4.170) 

y¡ (k) = T Lx, (x)] (4.174) 
where 

к) = Іх... (К) x Cole (4.18a) 
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a vector of 
yu CR? = Уха (К) 


a vector of 
subsystems. 


inputs to P linear subsystems; 


Y 
ge) CAE 


inputs to Q nonlinear memoryless subsystems; 
y, p(k) I! (4.180) 
БР і 
outputs from P linear subsystems; 
а Yugi? J Cie isd) 


outputs from Q nonlinear memoryless 
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Figure 4.8. A General Nonlinear System 
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The z - transforms of these signals are also defined as 

X; Cz), XX CZ), Y; (22, Y Ez). Г» Г.» C, and C, are matrices 
whose elements are either 0, -l or +l indicating the 
interconnections of the various subsystems. T[.] and F[.] 


are diagonal matrices specifing the linear and nonlinear 


subsystems as follows. 


Y; Cz) = T(z) Х (2) (4.19a) 
where A,(z) са 
i= C2) pe 
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This overall system given by Figure 4.8 is adequate to re- 
present a broad class of systems comprised of intercon- 
nections of linear and nonlinear subsystems including 
cascades, parallel connections and feedback systems. 
Equations (4.19) and (4.20) assume all the subsystems 
are single input single output noninteracting systems. If 
desired, the collection of linear subsystems given by 
T(z) can readily be put into the general multichannel ARMA 
E or Section III.F to allow each output to be a function 
Of past values of all outputs, and past and present values 
SI ELL inputs. 

Alternate representations of these linear and non- 
linear subsystems are also useful. Equations (4.19) are 


equivalent to 


Y, (2) = LA. (z) JX, Cz) + LB: (z) JY, (z) (uH 


En the time domain 


y, OO < La; CK) ]*x, Ck) + Lb: OO J*y, (к) (4.21b) 


Here * represents convolution and is carried out in the 
same sense as matrix multiplication and the matrices. 

La. (k) J, Lb: CK21, ГА; (2) and [B.(z)] are diagonal matrices 
whose i-th entries are the time domain functions а. (к) ог 
р; (©) or the polynomials А; (2) ор B.(z). The time domain 
manctions а. (К) апа b; (k) are the inverse transforms of the 


polynomials А. (г) апа В. (2). This can also be represented 
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in nonrecursive form as 
У (К) = ГВ; (К) 1х (х) i223) 


where Lh. Ck) ] is a diagonal matrix of impulse responses 


defined by 
И» h-(n) 6 (k=n) (4.22b) 
l n=0 1 


The nonlinear systems can also be represented in terms of 
inverse functions assuming they exist over the necessary 


ranges of the variables so that 


ES 5 
x (kK) EE Cy CK) J = 


егі 
F4 Ly O0] 
= ` (4.23) 


o 
Fo [Yugo 09 1 


So that equations (4.17) are iterative it is 
necessary that no delay free loops exist in the system of 
Figure 4.8, A necessary and sufficient condition for the 
absence of delay free loops is developed in Appendix H and 


requires that the terms of the determinant of the matrix 
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а = Бам Г, ltz) ii (4,24) 


Z> 


КО Шат its principal minors must be Zero. 
The various signal vectors and equations (4.17) can 


now be partitioned as 





4 
x; (kK) Па [уь [УСК С 
= + uCk) (+,25а) 
xy (K) ШЕП yx C 1b 
Da pop yp (kK) Ca] 
= + u(k) (4. 25b) 
Ге Га |У) сы 
p Баш НЕ Г. xy OÓ 
- ee 
оо ЕТО ЕТО ESOS 
уу (К) DNE an xf Ck) 
5 (4.254) 
ШОЛ ето та оо 


Equation (4.25c) can also be written in terms of inverse 
functions as in equation (4.23), and (4.25d) can be written 
using either recursive or nonrecursive representations of 

the linear systems as was done in (4,21) and (4,22). The 
single primed signal vectors in equations (4.25) represent 
those signals which are observed and the double primed signals 


are those which are not. It is assumed that all input signals 
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ШОО are Observed. It is possible to rewrite equations 
(4,25) as a single composite matrix equation as done in 
(4.263). As written, Е ово is an infinite memory or 
model order representation because of the presence of the 
IS operators but a finite memory version can be written 


by expressing rows 4 and 8 as 


уу (К) = а (К) ху (к) + Б (Ку (К) + а, (К)*ху (К) 

Nc S (4.266) 
у"(Кк) а (к)°х! (к) + ъЪ (к)у!(®) + а (к) х" (К) 

+ bà OO yr CK) (4.260) 


IN ira and seventh rows can also be written in terms of 


Bere functions if desired as 


es -1 -1 t 
КОО) Е УТ + Е. бу] (ц.26а) 
x(k) = Frifyt(e)] + FIALy"(k)] TED 
—N =c N -4 =} 


Now the problem of writing a system of iterative equations 
for the observed signals in terms of only observed signals 
consists of rewriting equation (4.26a) so that the upper 
ШЕПС Y by 5 partition is a null matrix. If this can be 


done, the general form of the nonlinear ARMA model can be 
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(х) н 0 Е 0 | UE NT 0 0 CD RA 

(х) 1А 0 0 0 sCOU | 0 б ет о] lox iA 

COD NX Oe ae 0 о | Da o ы 

(A) nX an 0 0 о | T 0 Oe SES Cine 

(892:) |-----| |-------------------- --------- ise Seen eee ee - |----- 
(х) {А 0 ОЯ о | От 0 0 (х) 1А 

(х) 1 0 0 0 (>) Ч | 0 o a 0 (х) 1 

OD TX одаи Ст 0 о | пб o 27. OR CD = 

= qT РТ eile OD Ix 
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used to identify the composite effects of the operators 
appearing in the upper left 5 by 5 partition. In some cases 
EES 111 опіу be possible of the infinite memory version of 
the model (with h(k)* operators) is used and in other cases 
a finite memory model will suffice, 

The process of determining whether infinite or finite 
memory nonlinear ARMA models can be used to identify a given 
system is illustrated for two examples in Appendix I. First 
a system consisting of a cascade of linear and nonlinear 
systems is considered. Then a model of the tracking behavior 


of a phase locked loop is put in the form of a nonlinear ARMA 


representation. 
2. Lattice Filter Methods for the Nonlinear ARMA Model 


As was the case with Volterra modeling, lattice 
filter methods can be applied to the nonlinear ARMA modeling 
problem if the regular form of the kernels is used. A 
family of signals is defined as nonlinear combinations of 


delayed values of y(k) and u(k) as follows. 


E. E ә (ІСІ) u(k-h, Ju(k-h,-h,)...u(k-h,-...-h,) 
НР) 
Yn, E = ж o, A) 
..-y(k-1-h,-...-h_) (%+.27Ь) 
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u y "one kh cuok-h,)-..u(k-h4-...-h. ) 
hj. A "hype hy, 2 2 р 


POLIS 
EI MC a (4.270) 
or alternately 
in, ch. “Uh КОО MEO анро -y(k=1-h,=...-h ) 
q q атр 2 2 а 
ОКЕ 350) 
(4.27d) 


With finite summations, the regular form of equation (4.14) 


becomes 
N N N 
P,P DB De 
5. > >, а Ch, n? un, 2 e 
0 h,=0 h,=0 
(4.284) 
Equation (4.15) becomes 
M M M 
q-q S 4, 1 
ее. bth, па) Прат 
eae h,=0 h,=0 
ЗБ) 
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and equations (4.16) become 


И L 
р+а,р+а рта, 2 p*q,i 
> > E Doc ES 22: Pis 
ШЕ 1-9 һ250 h, =0 
y (k-n) (4.280) 
ое о 1. 
9d 
Е i ii 
а+р,9+р а+р, 2 ODL 
D» Ср Ла ер) СЯ 
nt h,=0 h, =0 
u (Кер) (ц.28а) 
Ras" Paro ЗЇ 


These terms can be viewed as the summation of the outputs 

of a large number of linear filters whose inputs are the 
mals defined in equations (4.27). In the context of 
multichannel filtering, each of these filters can be con- 
sidered as a single channel and each of the input signals 
can be associated with one of the channels. Since the 
lattice models use the same amount of memory in each channel, 
the upper limits on all summations over hi will be made the 
ES N, 1 = DN = L *q,i = N.). The upper limits on the 
summationsover the other indicies determine the number of 


channels required. 
122 





For a quadratic nonlinear system, the present value 


of the system output is estimated as 


N, Nr Ny 
уб) = Denon) > 3 В 
DEI h,=0  h,=0 
“a Ый» ү 
5 Жо DD» > ео 
h, =1 h,=0 h,=0 
ШО 21 
+ У 2 6 E 5 266527 (4.29) 
А 2 1 2 
һ,50 h,=0 


where uy (К) = u(k)u(k-h,) ЗУ 
2 


B - у(к-І1)у(к-1-Һһ,) 
and u УҺ (к) = ulk)y(k-1-h,). Signal and coefficient vectors 
2 


can be defined for the various quantities in equation (4.29) 


following the conventions established earlier; for example 


T 
У(К) = Гу(К-1) ... УСК-М, ) 1 (ч.30а) 
Е [ y 
(С) = Ly, (k) y, (k-N,) (4.30b) 
E h, h, | 
апа 
i 
b = [b 01) ... b,CN,)] (4.300) 
Р Т 
Pn, * [5.(0,k,) ... bI, ,h5)] (ч.304) 
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Embedding all these vectors into single data and coefficient 


vectors, equation (4.29) becomes 


ую = X," da (4.31a) 
where 
i | 
он | ик у (ют 
1 У = Z9 : Ум, 
i ! 
EE iu ao 
| | 2 
i 
PIE | + T 
UY (KIT |... | A 1 SID) 
| | 
апа 
Т г. | Т 
Е a |b ; ...4 Da 
„+ | | M Т 
= | = 
0 | | Мо? 
T n 
| | 
e | | EH ] (4.310) 


The minimum mean square equation error solution for the 


model coefficient vector а, is given by 


2 d 
c£ X100 X400) 4, 7 et X4 GO yOO! (4.32) 
where the equation error is defined as 


e(k) = y(k) = y(k) 
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While equation (4.32) is similar in form to equations 


(A.7a) and (A.26a) for AR and MA models, 


lt lacks the nec- 


essary form for the application of the Levinson algorithm 


because of the fact that some component data and coefficient 


Meetors are N.-vectors while others are (N,+1)-vectors. 


1 


This is the same problem that arose in linear ARMA modeling 


and can be handled here in a similar manner. Defining 


p= [1 -a (0)  - b (0,0) 


- а (0,0) 
D 


- 710,0) 


and 


x(k) = Cy(k) uk) уш) 


4,40%) 


uyg lk) 


-b(0,M,,) 


-а (0,8,.) 


Т 


ze SC cr 


rl 22 


са) 


(к) 


u (k) 
22 


uy (k) J 
22 


(4733D) 


ШОС mines that the coefficients in y can be estimated in 


Some other manner, the MMSE solution for the remaining 


model coefficients becomes 
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ЕСХ (к) Хак) d = e{ X(k) x Got) ф (4.34) 


Here X (x) a E med es Х со and а. wich all 
superscripts"*" removed from their component vectors (indi- 
cating they are all indexed from 1 to N, and are all М. - 
vectors). Note that just as in the linear ARMA case, the 
ШОШО епз їп y essentially correspond to gains on their 
respective channels. Comparing equation (4.34) with (A.7a) 


for a multichannel autoregression it is clear that d can be 


obtained from the multichannel AR solution by 
а= ру (4.35) 


where the signals in x(k) comprise the channels in the auto- 
regression. Therefore, either the Levinson algorithm or the 
lattice methods can be used to solve for D and from it and 
Ши пон! едре of y, the nonlinear ARMA model coefficients can 
be obtained. 

By analogy with equation (3.27) it also follows that 
the nonlinear ARMA equation error is related to the AR pre- 


metion error vector by 


e(k) = р? elk) Е 


so that 


elek) =y Py (4. 36D) 
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Mere Pp is the AR prediction error covariance matrix. The 
coefficients in can therefore be set to minimize the mean 
square value of the nonlinear ARMA equation error in (4.36b) 


resulting in a solution given by 


E 


= É болуа) 


PN2 PNN А ET 
Mere N is the total number of channels in the model 


AS 1 + 1 + (М.2%1) + М, +1) + (101) 


and the pe are the elements of the prediction error co- 
variance matrix. It is readily apparent that the linear 
ARMA model and its solution via the Levinson algorithm or 
lattice methods are a special case of this formulation of 


the nonlinear ARMA model just as one would expect. 
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ЕЕ ВС В CONCELUSTONS; AND OPEN QUESTIONS 


In the previous four chapters, existing methods for AR 
and MA modeling were reviewed and from them new methods for 
linear and nonlinear ARMA modeling were developed. Here, 
two applications for these new methods in reduced order 
modeling and modeling for fault detection and diagnosis 
are examined briefly. Then the results of this research 
are summarized, conclusions drawn and significant open 
questions for the continuation and extension of this work 


EE listed. 


A. APPLICATIONS 
1. Reduced Order Modeling 

Oftentimes, complex physical systems, both linear 
and nonlinear, can be approximated quite closely using 
Simple models. The lattice solution methods developed here 
provide a very natural and efficient means of determining 
reduced order models for complex, high order systems especially 
in the case of linear systems. In Chapter III it was argued 
that for linear systems it is reasonable to build up ARMA 
models by simultaneously incrementing the order of the numera- 
tor and denominator polynomials as the lattice method does. 
When this method is used to build up a given order model, 
all lower order models and their mean square values of equation 


error are readily obtained as well (the only additional 





calculations needed are for the MSE and the gain term or 

ШЕП тшасгіх in the multichannel case) making it easy to 
compare the various models and decide if reduced order models 
provide sufficient accuracy. 

Consider for example the seventh order system whose 
characteristics are listed in Table 5.1. The magnitude 
Spectra of second, third, sixth and seventh order lattice 
models obtained using batch processing with 4000 point 
averages are compared to that of the system in Figure 5.1. 

It 1s apparent that a second order model is unable to 
approximate the system well, however a third order model does 
provide a good approximation. Furthermore, increasing the 
model order to four, five and six fails to significantly 
improve its performance as evidenced by the sixth order 

plot in Figure 5.1c. The model accuracy is not significantly 
increased until seventh order (which corresponds to the order 
of the actual system) when a very good fit is achieved. This 
is further illustrated by the plot of the normalized mean 
Square equation error as a function of model order shown in 
Figure 5.2. The cost drops rapidly going from second to 
third order but then fails to decrease significantly until 
seventh order is reached. The roots of the system and the 
various order models are also plotted in Figure 5.3. 

The benefits of the lattice method for reduced order 
honlinear ARMA modeling are not quite so pronounced however, 
Since adding extra stages to the lattice corresponds to 


allowing the kernels to grow in only one of their multiple 
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TABLE Sal 


SYSTIEN E 


Mee FUNCTICN COEFFICIENTS 
ПЕКСМ | КАТОК 


NUMERATOR 
A(C)* 1.C0000 
А(1): -2.16510 
А(2)= 1.5625С 


А( 2): 0.0 

А (4): 0.0 

А (5)= 0.0 

А (6)= С.0 

А (1) = 0.0 
ZEROS 

RE IM 


1.CE25C O». 6250C 
1.C5250 -0.62500 


с. С 0.0 
С. С 0.0 
С.О 0.0 
С.С 0.0 
С.С 0.0 


B(1)= 
B(2)= 
B(3)= 
e(4)= 
В(5) = 
в(6) = 
e(7)= 


ROAT LOCATTONS 


RE 
0.96000 
0.65282 
0.69282 
0:2941 
0422541 


-0.28750 


Ше) 


О 


2. 
-2. 


C. 


-C. 
C. 


18950 
21260 
51740 


=C.25250 


45144 
23462 
69321 


11} 
C^» 


IM 
0.0 
С.4009С 
-0.4000С 
0.557175 
20369518 
0.4979€ 
-0.497%6 
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Figure 5.2. Mean square value of equation error (as 
a percentage of the mean square value of 
system output) vs, model order for lattice 
models of system E obtained using batch 
processing and 4000 point averages. 
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dimensions. (The linear "kernel" had only one dimension.) 
Here some a priori information would be useful in determining 
the values to prespecify for the boundaries of the kernels 
in the other dimensions. Still however, when compared to a 
brute force matrix inversion approach where all kernel 
boundaries must be prespecified, the lattice method provides 
a viable alternative in obtaining reduced order models 
especially for low degree systems. 
2. Modeling For Fault Detection 
The problem of fault detection and possible diagnosis 
can be formulated as follows. 
a. Obtain a parameterization that describes the 
current functioning of the system under test. 
b. From this parameterization, determine if the 
S7secmetserumerloning normally or if a fault 
has occurred by comparison with a fault dictionary. 
It is the first part of this problem that has been addressed 
in this work. The parameterization can be as simple as 
sampled measurements of the response to specific inputs 
however, the large volume of data that would generally be 
involved in such an approach would greatly complicate the 
second part of the problem. A more efficient approach in 
terms of utilization of parameters is to model the system 
and use the model parameters as a description of its 
current functioning. 
For linear systems, ARMA models provide a very 


general framework with a number of possible parameter sets. 
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Three candidates are polynomial coefficients, root locations 
and lattice reflection coefficients with the latter offering 
many advantages. In addition to the advantages demonstrated 
by the experimental results of Chapter III, reflection 
coefficients provide a very effective and methodical way to 
build up knowledge of a system's characteristics. As model 
order is increased to more accurately represent a system, 
reflection coefficients already determined don't change 
making them ideal candidates for use in a dictionary lookup 
scheme (a characteristic not shared by the other parame- 
terizations). This is made more important since reduced 
order models could be adequate to detect and perhaps diagnose 
some faults, especially catastrophic ones. While more 
parameters are required when reflection coefficients are 
used, these same coefficients also provide all reduced order 
models. For a single channel ARMA example, 8N reflection 
coefficients and N gains provide all models from order 1 to 
N while N^*2N parameters would be required using either 
polynomial coefficients or roots to provide the same infor- 
mation.  (6N reflection coefficients are needed if the input 


EM te noise since K,.-k 00.5) 


1D 727 
A similar argument could be made for the use of 
MN Ice reflection coefficients with the nonlinear ARMA 


model for fault detection and diagnosis of nonlinear systems. 





E CONCLUSIONS AND OPEN QUESTIONS 

The purpose of this research was to extend existing 
theories and methods in the modeling of linear and nonlinear 
systems to broader, more general types of models. After a 
discussion of available results in AR and MA modeling of 
linear systems with particular emphasis on the Levinson 
algorithm and lattice filter methods, model transition 
formulas were developed to relate the more general ARMA 
model for linear systems to the AR and MA models. It was 
shown that with suitable assumptions, the ARMA model solu- 
tion could also be obtained recursively in order using 
either a modified Levinson algorithm or lattice filter 
methods. These results were developed extensively in both 
theory and practice for single channel linear ARMA modeling 
with experimental verification of both the batch processing 
and adaptive lattice methods presented. Portions of these 
results have already been published. [Refs. 65 and 66] 
The theory was also developed to generalize these results 
to the multichannel ARMA case. 

Based on the simulation results it was concluded that 
the lattice methods offer the following advantages over a 
conventional brute force matrix inversion approach to ARMA 
modeling using windowed correlation estimates. 

ime ron Short runs of data the batch lattice methods 

provide much more accurate results than the brute 


force method. 


1.98 





2. The batch lattice method performs much better than 
the brute force method when the system is over- 
modeled, 

3. The MSE as a function of model order is well behaved 
for the lattice method. 

^. The adaptive lattice method has difficulty tracking 
zero valued overmodeled parameters. 

The cost of these advantages is the extra computational 
burden of passing the data through the lattice filter during 
the modeling process. 

In the discussion of nonlinear system models the Volterra 
model was considered as a nonlinear extension of MA modeling 
and it was shown that lattice methods а be used to obtain 
the model solution if the problem was recast, using the 
regular form of the Volterra kernels. Then the new nonlinear 
ARMA model was considered and it was shown that this repre- 
sentation in some cases solves the problem of requiring a 
very large number of model parameters encountered in Volterra 
modeling. Then lattice methods were developed for the non- 
linear ARMA problem and it was shown that the linear ARMA 
techniques presented earlier are a speical case of the non- 
linear ARMA methods. For both types of nonlinear models, 
the recursive in order nature of the lattice methods was 
shown to allow the various model kernels to grow in one 
dimension while holding their boundaries fixed at pre- 
specified values in the other dimensions. The use of the 


model nonlinear ARMA was also illustrated with two examples 
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and a nonlinear ARMA model was proposed for the tracking 


behavior of a phase locked loop. 


several significant questions remain for the continuation 


w extension of this work and are listed here, 


Ir: 


Stability of the linear and nonlinear ARMA models 
must be considered. In the linear problem, stability 
1s dependent on the roots of the demoninator poly- 
nomial of the synthesis model. The methods 

developed do not guarantee stability of the resulting 
model. (This was not found to be a problem in 
practice, however, unless extremely short runs of 

data in the range of 30 to 50 samples were used. 

Even then, model instability was not a frequent 
occurrence.) Stability for the nonlinear ARMA model 
remains to be clearly defined. 

Input signal requirements in the nonlinear ARMA 
modeling process need to be investigated. In linear 
ARMA modeling the power spectrum of the input signal 
was found to play an important role. No requirements 
emerged however, on the probability density function 
(pdf) of the input. In nonlinear systems where the 
behavior is inherently level dependent, it is in- 
tuitively appealing to use an input signal with a 
flat power spectrum across the frequency range of 
interest and whose amplitude is uniformly distributed 


over the range of interest. 
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The inability of the adaptive lattice method to 

track zero valued overmodeled parameters is an 
interesting problem warranting further consideration. 
If some means of detecting the degeneration of the 
cost surface towards an infinite trough can be found, 
the problem could be remedied by simply reseting the 
appropriate parameters to zero. 

Experimental experience needs to be gained with the 
nonlinear ARMA model itself and the lattice methods 
developed for it. 

The characteritics of the lattice solution methods 
need to be further quantified to gain a comprehensive 
understanding of how and why it performs as it does, 
Also, further comparisons should be made between 

the lattice methods and conventional methods. Some 
comparisons were made here for batch processing 
methods but only a rectangular window was used on 
the data. Comparisons should be made using other 
types of window functions in the brute force method 
and the adaptive lattice method should also be com- 
pared with a conventional LMS adaptive algorithm 
applied to the equation error model in which all of 
the a(i) and b(i) coefficients are adapted simul- 
taneously. 

In the adaptive lattice method, scaling of the lattice 
input signals needs to be investigated. It was 


noted that the ratio of largest to smallest eigenvalue 


zur) 





of the input autocorrelation matrix was related to 
the speed of convergence of the adaptive algorithm. 


For the first lattice stage this matrix is 


|o? (o M ES 


mer NEU) as 


If the system has a high gain such that the mean 
square value of the output y(k) is very much greater 
than that of the input u(k), convergence will be 
slow, This could be remedied by implementing an 
adaptive scaling scheme at the lattice input (perhaps 
similar to the first order low pass filter estimates 


used for adaptive step size normalization). 
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ACPENDIX A 


Alternate Multichannel Model Forms 


Multichannel generalization of the MA and AR models were 
discussed in Chapter II along with their solution via the 
Levinson algorithm. Here the multichannel models and the 
Levinson algorithm are developed in an alternate form more 
compatible with other linear and nonlinear modeling problems 
to be explored later. 

Consider first an N-th order, Qg-channel, AR model where 


Meemenrrent value of the signal vector x(k) 


- 


a 


[x Ck) Mec С) | AID 


E 
x(k) 


T 


LX. Cz) БЕ ун (А.1Ь) 


Х( 2) 


mS be predicted from a weighted combination of N past 
values of each of the component signals. For each signal 


this estimate can be written as 


So М 
> З ОКП) САК 2а) 
= Ра 1) 1. 


СК) = 
J НЕ! 
or 
Чо 
Y. Ж а..(2) Х. (2) 22 
Хх. (к) E ij 2 1 2 
where N 
= =]] 
БЕС? = 2 2 (А.2с) 
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Define an N-vector for each of the 99 channels to contain 
their required time histories as 


ДЕ 
ии е кейк 1) 7, жк. (Ck-N)] (доза) 
= gi т 


ШОО < 1 < Oo and embed all of these vectors into a single 
NQ,-data vector given by 


T 


` T T 
AK) = [x (k) wt xe URP (A.3b) 
£z» = “90 


Define a М9, x Qg matris of weights as 


D E 4,4 a “19, (A.4a) 


2Q 1 9%) 


where the N-vectors ШЕ are given by 


Ж JT. 
age = Ld: (1) E dj (N) J (A.4b) 

and contain the coefficients of the polynomials SE 

These polynomials can be combined in a matrix polynomial 


defined by 


(A.4c) 


Dez) 


— 
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ши и спеве definitions, the N-th order prediction of x(k) 


and its associated prediction error vector becomes 

x(k)? = Xx)" D | (А.5а) 

e(k)? = x(k)? = x(k)! (А.5Ъ) 
or in the transform domain 

E(z) = x(z)" [I - D(z] (А.5с) 
Comparison of equation (A.5c) and (2.44c) show that this 
matrix polynomial differs from the more generally used form 
of (2.44c) by a transposition. The coefficient matrix D 
can be found by minimizing the trace of the prediction 
ево Covariance matrix 


Р = є{е(К) е(к) 1) (A.6) 


leading to a system of linear equations given by 


«Хо Хо D = eX xo (A.7a) 
or 
E. | 3, xo = Tx x | 
n % 
р = : 
“Xo = u В 0 хх T. x 
>. 9 0 79 | Qo 1 So 
(A. 7b) 
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Adopting a shorthand notation, this becomes 


Ер = г (A. 7c) 
As in the case of a single channel autoregression, the 
multichannel generalization of the Levinson algorithm also 
requires that the backward prediction problem of estimating 
X(k-N) backward in time from the values of x(k-N*1) through 
E SS be solved simultaneously. Using an overscore to 
indicate quantities associated with the backward problem it 


hollows that 


== T ОТ 
Хю D = x(k) (09935 
T HN ST 
е(к) = x(k-N) - x(k) (A.8b) 
DI 
p x = 
ва) = XCz) 2 Nez - Diz) (А.8с) 
== ат " 8 
where ХКК) = (х, (К) e XQ 0071 (A.8d) 
0 
m 
Е СЕМ)... х. (к) (A.8e) 
FTL al 1 
- г E 
and Р = іа, dio. (A.8£) 
а. 2,23 
ELE р 
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a) [di (2) es 42080. (A. 8g) 
with Diz) = " a (A,8h) 
d ZI X 
EU 229% 
N 
and d..(Z) = MA (A.81) 
23 АЕ т) 


Setting the coefficients of this backward predictor to mini- 
mize the trace of the backward prediction error covariance 
matrix 


P = efe(k) S(x).) (А.9) 


beads to a MMSE solution given by 


ET ете ls 
e(Xao Xao ) D - ={Х како!) (A.10a) 
. = = T ۳ 
and since Е {х, (К) X (x) } S 
and € (x; (Kk) x; (k-N)} = 


mars can be written as 
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Em > 





үз 
— Al ЕТО а E 
0 0 
D - 
RC X "ui RC x Ех x I x 
Qg 1 99 90 E Чо 99 
(A.10b) 


Adopting a shorthand notation, this becomes 


R 


I2 
“1 


(A EL) 


At this point it is important to take note of a subtle 
difference between the single and multichannel problems that 
has arisen. While the single channel equivalent of equation 
(A.10b) for the backward predictor was not written previously, 


ШЕ 15 Obviously 


R E d E up (A оа) 
Se] 


End since 


R = К (A.12b) 


it follows that the single channel forward and backward 
prediction problems have exactly the same solution (as was 
found to be the case in Section II.D. when the AR prediction 


error lattice was developed). This fact is responsible for 
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a number of simplifications in the single channel case: 

a.) Only a single reflection coefficient ge E Was 
needed to link the n-th and (ntl)-st order solu- 
tions for both the backward and forward problems 
КЕП Шота (2255) апа (2.,25с)). 

ЕР) |58(22)| = |В(2)| 

с.) є{е(к)°} = є{е(К)^} 

d.) The development of the Burg method and the Itakaura- 
Saluoumetnod for Calculating the reflection 
coefficients are a direct consequence of c above. 

Unfortunately however, none of these simplifications carry 

over into the more general multichannel AR problem because 

by comparing equations (A.7b) and (A.10b) it is evident 

that in general D 2 D 
Peoceecing aS in the single channel case, it is shown 

in Appendix B that because of the structure in the correlation 


matricies R and R , equations (A.7c) and (A.11) can be 


re-written as 


В =P (A.13a) 


and 


R F = Ша (ДО 5) 


where ЕЕ. p and е are obtained from D, D, I and 


I respectively by taking their component vectors and turning 


them upside down in place; i.e. 
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o 
o 








Een Bus 
Е = ; ` (А. 1ца) 
f fT 
ELE 0% 
where 
m 
а >] (A.14b) 
253 nij ij 
and 
P- E Р (A. luc) 
ЕЕЕ m E 
Б, X C X 
| 95 1 05 90 
with 
ар Ben (A.14d) 
i XiX Eu x. 
er 2 


As will soon be evident, the relationships of equations (A.13) 
form the cornerstone for the Levinson algorithm and are made 


possible because the blocks comprising the R and R matricies 


are themselves Toeplitz and because R = e í 
= ER 


Zssume that the (n*l)-st order solutions can be related 


Шле n-th order solutions by 


400) Боз 
=i} =i 5 
gor. TE, Аве 
0 (ntl) 
| a 
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and 





к x = 
BE Sn ---- (А.15Ь) 
J | 0 [porn 
J a 
Embedding the vectors ei? and E and the coefficients 
„(nt1) and к( 71) into matrices designated c and cim 
ij 1 5 = = 


к(п+1) and gem 


Laws cH пори от them in the (n l)-st order 


problem it is shown in Appendix C that 


— (n) 


К ы” (A.16a) 
= و‎ Ee ща (A.16b) 
(n*1) EO D B Din)! Hin) 
NM Ие Dei пк QGrn-Z "D" 


(А. 16с) 


T -1 ЛЕ T— 
(n+1) ER ET BD СВ. (п+1) E Шел 


о 


LAI 
T 


(А.164) 


Also the inverted terms on the right hand side of equations 
(A.16c) and (A.16d) are just the backward and forward pre- 
diction error covariance matrices for the optimum n-th 


order models so that (A.16c) and (A.16d) become 


21l1 


E S S S 





Di — Ще 
«(n+*1) = E [В (n+1)- 2° DB CID 


= Ion) Hi) 
[R „(n+) -2° DI (A.17b) 


mont) 
As in the case of the single channel problem, the prediction 


error covariance matrices obey the recursions. 


E E il), (A.18a) 


pint) p ome nS | (A.18b) 


ято 


completing the multichannel и of the Levinson 
algorithm for AR models. 

Bemparing equations (A.16), (A.17) and (4.18) with their 
mete channel counterparts (2.18) and (2.19) it is clear 
that the multichannel Levinson algorithm simply represents 
a matrix algebra generalization of the single channel al- 
ШІСІ. Once again, predictors of all orders 0 < n < N are 
obtained in the process of finding the N-th order predictor 
along with all their prediction error covariance matrices 
and the overall minimization requiring the inversion of 
a QA x QgN matrix is replaced by a sequence of N minimiza- 
tions, each requireing the inversion of two Qg Ж 00 matrices 


Да 








Using the relationships given in equations (A.15) and 
(A.16), successive orders of matrix polynomials can also 


be related to one another by 


EBD)? 0951-z- 527" cr- DO? cz; t» 


— 


(A.19a) 


BE [up (2) 1527 Та 0 1- ОС) CD к 


(A.19b) 


Premultiplying both sides of these equations by X(z), trans- 
posing, and transforming into the time domain provides rela- 


tionships between the prediction error signals at each stage. 


(п) 


ВШ = ето 


iU 
- ru e) 1) (A.20a) 


Ш) Gp ceni oS 
e DURS e 


ПОС е9, е` '(k) (A.20b) 


Recognizing that for a zeroth order prediction, the forward 
and backward prediction error vectors are just the input 
vector itself, it follows that 

(0) 


е Ga = 6 


BO Sex (A.20c) 


and the multichannel equivalents of equations (2.28) have 


been obtained. 
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ШО Con IGE che N-th order multichannel MA model in 
which the current value of a Qg-vector Of output signals 
етнос to be predicted from the present value and N past 


values of a Q,-vector гори signals x(k). 


0. М 
у. (К) = 9 9 d..(n) x, (k-n) m 
J T) 1. 
о 
Qi 
ве д b 
bez) = dq cz) (A.21b) 
J 1] 1 
i=l 
where N 
at. (z) = 2 d..(n)z P А ще) 
ЩЕ) 4,2) 
n=0 


Using a superscript "+" to indicate the fact that the vectors 
are indexed from 0 to N rather than from 1 to N, define 
(N+1)-vectors for each of the input channels to contain their 
required time histories as 


T 
x: (X) = Lx. CK) Bt х.(К-п)] (А.22а) 


and embed these vectors into а single Q,(N+1)-data vector 


given by 


X'ao - [x1 O07 Ш ох 011 (A.22b) 
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Define а Q.(N+1) X 9% matrix of weights as 


2. + + 
D = 9-1 м 310, $539 
+ + 
а а 
p 741401 


where the (N+1)-vectors di are given by 


MEC а орт (A.23b) 
55 1) 1) 
and contain the coefficients of the polynomial dii C2). These 


polynomials can be combined into a single matrix polynomial 


given by 
+ + + 
О => = jt» РР ч Сое) 
= + 
d (2) eas d m) 
| ат en 


cese definitions, the N-th order МА prediction of y(k) 


is given by 


yGol « X'aoo? D' (A. 24a) 
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or in the transform domain, 

Yo) = xc»! D'co (A. 214b) 
Defining a prediction error vector as 

eg 007 * yao! - yao! (A.25a) 


: ere : т К s 
ШІ сесііпр the coefficients in D So minimize the trace 


of the prediction error covariance matrix 


_ T 
Р = бе (Юе (КО) (А.25Ъ) 


results in a solution given by 


є{Х ao X'aob D' - etX'oo yao! (A. 26a) 
or 
+ Se ay ке 
X1 X1 Е X1Y1 т 
р” - я 
Е + xd у Е + ig = + y = + 
5% E 01 91 Qi i 0:7% 


(A.26b) 
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Adopting a shorthand notation this becomes 
+ 
Ер: т" (А. 26с) 


Assume a relationship between the components of the 


(n+1)-st and n-th order solutions given by 


515 415 
aU И 7 Ce 
0 A | 
1) 
Embedding the vectors E and the coefficients RI 


into matrices designated nsn? and а 


, and solving 
for them in the (n+1)-st order MA problem it is shown in 
Appendix D that 


jn*1 . p(n*D gt (A.28a) 


ES > T 
g(n+1) „ g(n+1) Ee Dr; ES 


—(n+1) —(ntl) 
where F 2 amd ри 


are matrices, that emerge 

from the backward prediction problem in a multichannel (n+1)-st 
Erler autoregression on the input signal vector x(k). Again, 
ШЕ 15 Clear that the multichannel MA solutions given by 


equations (A.28) are a matrix algebra generalization of 


equations (2.23) for the single channel MA model. 
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То relate successive order MA matrix polynomials to one 


another, equations (A,27) and (A.28a) can be used to write 
ШИСИ ST Dr Dayo (A.29) 


where the second term on the right hand side that premul- 
tiplies G is the backward prediction error matrix polynomial 
ШЕШ “Пе autoregression on the input signal x(k). Pre- 
multiplying by ХЕ (г), transposing, and transforming into 

the time domain results in the multichannel equivalent of 


 Шастоп (2.37) 


An) Со, "gan 


a nU = УП) +6 


(к) (A.30) 


and completes the derivation of the recursive in order 


Pomurtions for the multichannel MA model. 
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APPENDIX B 


A Key Relationship For The Multichannel Levinson Algorithm 


Equation (A.13) is a key relationship in the development 
of the Levinson algorithm responsible for much of the 
algorithm's simplicity. Without this relationship, more 
frame just the K and K matrices would be needed to obtain 
the new model from the previous model. 

In equation (A.7), consider the multiplication of the 
i-th block row of R by the j-th block column of D To 
form n 


i j 


R С К A ST (Bes) 
„X; X 1] ES E 


In particular, consider a general term on the left hand 


side Of equation (B.1) in detail. 


i (0) m ES (1-М) E К „61 

E im om M сле т 
I (N-1) R (0) а CN) R P (N) 

м л mJ = 
Б 


Define upside down versions of the E and РР vectors 


ij 
a 
1) 


Т Е (Bia) 
eu] ° 


aS 


Съ) 
1) 


21.3 





о 
ii 


RTT : CB. 3b) 


Using these permuted vectors in equation (B.1) in place of 
ELE nd © vectors, the relationship is still satisfied if 
BEEN -r»ces are permuted as well. In particular, from 


Meme) it is evident that 


Re (0) De (N-1) E (N) 
m m 
m 
R y (1-N) к. ШІ) ТШ 
т Lm lj 
R QD 
ON S 
1 3 
о. = (B.4a) 
R (13 
KAR: 
T) 


or equivalently 


T T 
R IN EER pcs Е (В.Б) 
к.к. =13 Е — 9 XiX; 
Embedding the Sr and р. X. vectors into matrices designated 


Lg — 
F and // respectively and using the definition of R nu 


follows that 


ВЕ A E) 





Defining E and 2, „ as upside down versions of their 
| =< = 
corresponding vectors in D and I in equation (A.11) and 


e с) 


embedding them in matrices designated F and 0 it also 


follows from a similar development that 


КЕ. 


ISl 


ES) 


Equations (B.5) and (B.6) are essential to the develop- 
ment of the Levinson algorithm and are made possible by the 
ВЕРЕ structure of the block components of the R ana 


— 


R Matrices. 
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A PENDIX C 
The Multichannel Levinson Algorithm Solution 
For AR Modeling 
In the (n*1)-st order versions of equations (A.7) and 
(A.10), the component matrices that make up R and R een 


Bemuritten as follows 


| жак (-М)| 
(п) 
R | 
ent) 2x. X. 
x. EE | а. 
1 7) ИСЕ и іы 72.0 
Т 
СЕТТІ х о К, EOD 
i 15: \ 
| 
Е ри 
= ij | + J (oa) 
ee ss +---------- 
i | 
(п) 
E. | К Sc» 
jn 
and - 
T 
Nr) (nti)? в) | а. 
SUR = d |] T (C2) 
E Ox. EX. l 
и IN سسس‎ Ec 
Ku | 
= | St UD) 
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Additionally, 


p intl) = ml d ——— 


EIS AX 
25773 В 


СЕЗ) 


With these matrices and vectors written in this partitioned 
form, and with the relationships assumed in equations (A.15) 
between the n-th and (ntl)-st order solutions, the (n+1)-st 


order modeling equations become: 


Forward Model: 
R” D т Ru и E ОО аг: (C.4a) 


— T —, АТ 
¿DS Do | E» a " в (ok tb 4 к, (ntl) 


(C.4b) 


Backward Model: 


RO po , go» Е), ом кар . то 


. 
pd س‎ „р 


(5/52) 
T — urs c 


C55) 
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Equations (C.4a) and (C.5a) contain the n-th order modeling 
equations within them however, and therefore can be written 


as 


RS? eu a Jo ¿A (C.6a) 


p a " po Un S (C.6b) 


and applying the relationship developed in Appendix B 


(equations (B.5) and (B.6)) yields 


co)? TR E x On*1) (C.7a) 


В ко 


АІ 


КСО) 


us the = and € vectors have been found in terms of the 
known quanties f and f and the unknown K and K matrices. 
ВБ ЕцЕ11> equations (C.7) into (C.4b) and (C.5b) completes 


the solution resulting in expressions for the K and K 


matrices given by 


Eu е == 20% 
mtb 2 [R (0) - DS Dj [R,. | (n$2) 2? DE 


(Сева) 
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T -1 Т--- 
Kan = гк, (00 r ID. Св. „Сао? 9C? Ds 


ЕЕ) 


The terms on the right hand side of these equations that 
are inverted are just the backward and forward prediction 
error covariance matrices. for the optimum n-th order models 
given by 
(m) (п)! (п) 
г р (C.9a) 


хх 


| o 
ii 


L 
а 
— 

il 


Bey Seep. 
R (o) - rm? pt (ono 


Using equations (C.3), (A.15) and (C.7), the forward pre- 
diction error covariance matrix for the optimum (n*1)-st 


order model can therefore be written as 
eet) (m) m GO, e GO E GO, (n*1) 
El quc Mor к у 


Ban 


оо) (Come) 
—хх 


Та. 
Е Е спо k (С.10Ь) 


==“ 


dE == 
я Б оо КК 2 DN и 


02 с) 














pD. pM pr pD KOD] RR 


— 


and following a similar development for the backward pre- 


Aeon error covariance matrix results іп 


(n+1) 


Der _ tl), (C.11) 


P 


Ше 
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SPTENDIX D 
The Multichannel Levinson Algorithm Solution 
FO IA Modeling 
First note that because of the definitions of the x; (X) 
and x; (k) vectors (indexed from 0 to n and 1 to n respec- 
tively) in the n-th order models, the matrices БА үп 


the n-th order MA problem could also be written LJ in 


terms of x, (kK) and а 5) as (assuming stationarity) 


(и) _ ШІ) 
Ев = Ren. (Dora) 
X. fie 
1 3 
so that 
(п) 
т (n+1) 
R R (D.1b) 
ponto) ашт 
The components of апа Г in the (ntl)-st 


order MA modeling problem can be written in partitioned form 


as 
m. | —(п+1) 
к Be 
(n+1) er ini 
R - 
— + + T"—————————— ست م‎ (D. 2a) 
Ex. 
ае T 
RI 
Py x. О) 
ШЕНІ | 


2. 





and 


(n) 
r 
(п+1) Xy. 
n = 171 ЧБ DD 
EE EM 2250 
R o 


Using these partitioned forms, and the relationship in 
equation (A.27) between the n-th and (n*l)-st order solutions, 


the (ntl)=st order modeling equations become 


en) i) 
R+ D* 


(п) — (n) 
: Rt gen à осоЕ y + 


(Шла) 


n) 


par, a) an (п+1) 


+ R (0G 


а oul ТЕРЕ!) 


Eomacıon (D.32a) contains within it, the modeling equation for 


the n-th order MA model and using equation (D.1b) can be 


rewritten as 


Roth 0 +1) 2-2 отан Ds 


A comparison of this result with equation (C.6a) shows that 


| (n+2) M F (ntl), (ntl) (D.5) 
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F (atl 


where is the permuted version of the backward 


pemueion in an (ntl)-st order AR model of the input signal 


вое x(k). Furthermore, substituting equation (D.5) into 


(n*1) 


E c) results in a solution for G given by 


= T Зи 
т) СВ. (0) u m) E y 


Е T 
[Rynt pi] (0.6) 


Since 


a) got В ток) рсн) 


it follows that the inverted term on the right hand side 
of equation (D.6) is just the backward prediction error 
covariance matrix for the optimum (n+l)-st order AR model 


On the input signal vector x(k). 
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APPENDIX E 


Prony's Method For ARMA Modeling 


Prony's method [Refs. 8, 52 and 56] obtains a zero pole 
model for a system by matching the impulse responses of the 
System and model over the first M+N+l sample intervals where 
M and N are the orders of the model transfer function num- 
erator and denominator polynomials. Assume that a signal 
y(k) is available that represents the impulse response of 
a causal system and that a rational transfer function model 
for this systems is desired. Using a "^" to denote the 
model output and u(k) to denote the input to both the system 


and model, the model transfer function is given by 


И -п 
х DEO 
“СИЛ 21-0 = 
Hz) = Е = Bela, - BE) 
пы аля 
[| = 


For a unit sample input it follows that 
Bea Z2) AZ) (E) 


Equating like powers of z in this relationship results in 


N on) E Шоп М 
E. b(i) y(n-i) - So 
1=0 


0 else 


where b(0)=1. 


23:0 





Equating y(k) and y(k) over the interval 0 « k « M*N produces 


a set of M+N+1 equations which can be expressed in matrix 


Ferm ас 
[co | о 0 ЕП? ---- 
| b(1) 
СІ) | y (0) 0 LB ; 
Bm) 
E | у(М-1) у(М-2) B 
E E ас 
yon ІІ 001 
| 
| 
Қасы» y (D, 


or adopting a shorthand notation 


Г + m 
Y» Y, 5 0 
L B 


Assuming that the NxN matrix Y, 


mm ср -1 
а 15 y 


| 
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ad) 


а(М) 
0 
(Т д) 
ца 
(ESUD) 


Sino sineular 11 follows that 


VES Sa) 


ESD) 





The original Prony method goes on to form a partial 

fraction expansion and inverse oH on the model 
transfer function H(z) resulting in a model for the impulse 
response of the system given as a sum of complex exponentials. 
This in unnecessary here however, since a rational transfer 
function model is the form sought. 

Prony!s method inherently assumes that matching а 
sufficient portion of the impulse response of the system 
results in an accurate model but this is not necessarily 
the case unless the impulse response damps out quickly or 
unless the system can be represented exactly by a low order 
model. Otherwise a very high order model may be required 
to obtain sufficient accuracy. Other difficulties asso- 
ciated with this technique are: 

1.) The system impulse response must be available; 

2.) There are no built in mechanisms to test for or 

ensure stability of the model; 

3.) There is no averaging of noise in the data; 

4.) Only a small portion of the available data points 

(M+N+1) are actually used. 

These difficulties can be partially overcome by modifying 
the procedure to obtain an approximate match of the system 
and model responses over their entire duration rather than 
асс шассп оуег the first M+N+1 points. Consider 
equations (E.4) modified to include the entire signal y(k) 


For 0 CE >. 


372 








yo) 


y 
p I) 


г + | 
Xi Yi 
Ша Хз 


but with y. >; У. 


e@mation (E.7b) 


equations 


[ед 


о 


and 0 having an 


| 


[o 


(Е.ба) 


(Eo 


СЕЙТЕН 


CEB) 


infinite number of rows, 


will in general have no solution. In practice 





only a finite portion of the system impulse response y(k) 


can be considered but equation (E.7b) will still be incon- 


а Іп general. А least squares estimate of b can how- 


es T 
Eu be Obtained by minimizing e e where 


(E.8a) 


СЕВ) 


: Қ қ і + 
AN turn can be used in equation (E.7a) to finda . 


This approximate version of Prony's method is the one most 


commonly applied. 





A FPENDIX E 


Parabolic Surfaces In n-Dimensions 


Multi-dimensional parabolic surfaces are described by 


an equation of the form 


KA. 2K Dre (F.1) 


where: 


is the independent variable; 


к 


is a vector of dependent variables; 


|ж 


| > 


is a symmetric constant matrix; 


ІС 


1s a constant vector; 


Eo constant. 


О 


(х, b and c can also be considered as matricies with the 


posce of the right hand side set equal to the independent 


variable but the problem remains essentially unchanged.) 


Sempreting the square for nonsingular A this becomes 


= шта 


(БА?) 


= 
it 
m 
| 
I» 
Іс 
|» 
| 
| 
I 
| 
ке 
іс 
O 
| 
O” 
| > 
Іс 


EN ust fOr positive definite A it is clear that the minimum 


Кете ог y is obtained for x = A 


с - b A tp. It 1s also clear that nonzero values of b and c 


b and this minimum is 
simply raise and lower the surface and move the minimum 


Point away from x = 0. The shape of the surface (its 


relative concavity or flatness) is determined by the matrix 
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БЕ спе Guddratic term Of equation (F.1). Therefore, to 
study the shape of this surface, consider the simpler эго- 


blem with Бл сеш Zero - 
о (E35) 


One way to examine the relative flatness or concavity is to 
look at the locus of points on the surface for constant 
pues of y; in particular when yl. Recognize that A can 
Bde rewritten as QA ol where A is a diagonal matrix of 
eigenvalues and Q is a matrix whose columns are the unit 


ШЛЕ eigenvectors of A. Now (F.3) becomes 


[x 
КО 
| > 
| 


xc (F.4) 


and introducing a new set of variables w=Q x (which are 
Bey So rotation of the variables in x), equation (F.4) 


reduces to 
2 
Ac vus I = 1 (F.5) 
mE 


This equation describes an ellipsoid in n dimensions whose 
axes half lengths are given by 1//à; for 1 <i <n. This 
Pollows from letting all but one of the w's equal to zero 

and solving for the nonzero variable so that one point on the 


surface is for example на 7 ТИ УЛ: with Е Ши 


1 
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Mies Point 1S just a multiple of the first eigenvector of 
A so that in general, the axes of the ellipsoid point in 
the direction of the eigenvectors of A with half lengths 
given by the recriprocal of the square root of the eigen- 


values. 
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oe LG 


Multichannel ARMA Modeling 


Consider a system with Qg output signals Ly4 OO . « Yo (k)] 
Ü 


and Q; input signals ч, (к) E. Uo (k)]. The multichannel 
1 

ARMA analysis model forms an estimate of the present value 

of each output as a weighted combination of past values of 


all output signals and past and present values of all input 


signals. 
% N 
х u dy». ук: 
Y. (k) بے‎ Ts Ya i) 
Qi м | 
р » > a. (1) pcc (G.1) 
o : 


Define data vectors for all the input and output channels as 


_ т 
y,GO * [y,(k-1) ... y, (k-N)] (ера 
Mme) = fu (k) ... u (k-n)}* (G. 2b) 
ЕП п n 


and embed them into QgN and 0; (+1) vectors given by 


T mu 
Y = Ly, (x) 15: ы” 1 (G.3a) 
+ DES ТЕ 
U = lud щ © (G. 3b) 
e 





Define NQg x Qo and (N*1)0; x Qo matrices of coefficients 


given by 
- 
а Ва 
Чо 
B> ШЕР ша) 
b b. 
EM 0% 
and 
Г + + 
211 2m 
90 
At = |: (G. ub) 
E E 
Ril Bi 
where 
_ Т 
bs Db) к (уут (в.е) 
+ 
ai; = [a,,00) а) (в.а) 


With these definitions, the multichannel ARMA estimate for 


the vector of output signals becomes 


yao! = Cy) | и” 1 E (5.5) 


Forming a prediction error vector as 


eg) * yOO - yOO (G.6) 
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and setting.the model coefficients to minimize the trace of 
the prediction error covariance matrix results in a solu- 


tion given by 


= 
R R B R 
—YY mou E =Yy 
= (G.7) 
+ 
t à R 
ш ү UU U У | 


In the transform domain, the prediction error model is re- 


presented by 
E = Xtz)" [I-B(z)]-Ulz)" Alz) (G.8) 


where E, Y and U are the transforms of the error, output 


=a?) 
and input vectors and the coefficients of the i,j-th elements 
ENErhnecmatrix polynomials B(z) and A(z) are the elements of 


the vectors ЕВЕ апа а... The multichannel ARMA synthesis 


model is then given by 
Y'(z) = U(z)" A(z) [I-B(z) 7+ (С 9) 


with the matrix polynomial fraction serving as the generali- 


zation of the zero pole transfer function. 
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ОРЕХ Н 


Delay Free Loops 


To develop equation (4.24) guaranteeing the absence of 
delay free loops in the system of Figure 4.8, consider the 


equation for yy Oo. 
уу(К) = Ех (К) ] CER 


Since F[.] is a memoryless nonlinear function, proving that 
ху к) Ш5 пос а function of yu (X) is equivalent to proving 
that yu OO is not a function of itself, and therefore no 
delay free loops exist. From equations (4.17b) and (4.17d) 
with u(k)=0 it follows that 

И-П 


NOT = (H.2) 


2 ah 


and the coefficient of Yy at time k on the right hand side 
is given by 

a = Lim T, г) Г (Нез) 

Zo 

A nonzero ШЕ indicates a dependence of ху СК) upon Ying Se 
and clearly therefore all the main diagonal elements of a 
I be 2eno to avoid delay free loops. These elements are 
the terms of the (Q-1)-st principal minors of a, While this 
is a necessary condition it is not sufficient to avoid delay 


free loops since loops may exist through two or more signals. 
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The condition that EE Zr Ое x.) such that 
ES; < Ч апа 273), ensures that no delay free loop exist 
through 2 signals and is equivalent to requiring all terms 
ЕПС” (0-2)-па ргіпсіра! minors of a are zero. A term of 

a determinant [Ref. 63] is defined as the product of elements 
of the matrix taken one from each row and one from each 
column and the determinent of the matrix is the sum of all 
possible terms. It is clear therefore that every term 

must contain a cycle such as ES eur and must there- 
fore be zero. Since the determinant consists of every 
possible term, requiring that all terms of the determinants 
of the (Q-i)-th principal minors are zero ensures that no 
delay free loops exist through any combination of i signals. 
Examining all terms of the determinant of a and all its = 
principal minors ensures that all possible loops through 

the Q signals in Xv Ok) are examined. If any delay free loop 


exists, then at least one of the terms of one of the deter- 


Banants will be nonzero. 
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APPENDIX 1 


Nonlinear ARMA Modeling Examples 


The determination of memory requirements for the non- 
linear ARMA model as well as its applicability to systems 
consisting of interconnections of linear and memoryless non- 
linear subsystems is illustrated here using two examples. 
First a cascade of linear and nonlinear subsystems is con- 
Sidered. Then a real world example 1s considered and a 
nonlinear ARMA model is proposed for the tracking behavior 


of a phase locked loop. 


А. A CASCADE OF LINEAR AND NONLINEAR SUBSYSTEMS | 
Consider the system shown in Figure I.1 where the signals 
u(k) and Yr gk) are observed. In terms of the topology of 
Р; | ; . та: 
igure 4.8, seven signals can be identified On a X,» Ута, 


Уго» Хүү» Уу and u) however for convenience three of the 


seven equations in (4.26) which specify 


Хр СК) = u(k) Ci a) 
Xp (K) = yy1 (СК) Gl) 
Хуа O90 - y; 109 СА? 


will not be explicitly written. In this case, equation 


(4,26) becomes 
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Uk) 


ar > a -> «а» 





MZ) 


where the finite memory representations of T, (2) and T, (z) 
have been used. The objective now is to eliminate a, (k)* 
from the upper right partition. Using the equations of rows 
three and four, the first row can be rewritten as 


pes) =D, (K)* y, ,(k) +a, (k) *F, La, (k) uk) tb, (k)*y, 1 (O1 


Y12 
GE 


Estronally it is difficult to write equation (I.3) in the 
Operator matrix form of equation (I.2) because of the non- 
linear function however, the following representation is 
adopted. 


uy 





(h' I) 


© 
| | 
E 
ha 


© 
^ 
52 
— 
Q 


Г ab Ge ep ma Gee qe feet Ame Wb eee Өте» ue 


e 
e 


UE CUM EE 


т. «эш» < = = = سے‎ > = < = <= = > = > > > > <> > > > > <> > 


Т 0 
ыы Сы 
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No further reductions are possible to make the upper right 
Paeettlon a null matrix leading to the conclusion that for 

a finite memory nonlinear ARMA model to be appropriate, 
either b- (k) must be zero (the first linear system must be 
nonrecursive) or Y 1 (К) must also be observed. Alternately, 
ТЕ bi (X) is nonzero, an infinite memory representation can 
be used for the first linear system by replacing aj CC.) 
with h4 CO*(C.) and b, OFC) with zero in equation (I.H) 
indicating that an infinite memory nonlinear ARMA model is 


appropriate when only u(k) and Y1 2K} are observed. 


ОО NONLINEAR ARMA MODEL FOR A PHASE LOCKED LOOP 
A continuous time model for the tracking behavior of a 
phase locked loop [Ref. 55] is shown in Figure I.2 where: 
94 Ct) is the phase of the incoming signal 
9, Ct) is the estimate of the phase of the incoming signal 
e(t) is the phase error signal 
ШЕ) 15 the transfer function of the loop filter 


K and Ко are constants 





Figure I.2. A nonlinear model for the tracking 
behavior of a phase-locked loop. 
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NN cdel Ts nonlinear because of the sin function in the 
loop. Often the assumption is made that е(+) <<п/2 so 
Sas in e(t) 2 eCt) in which case a linearized model is 


obtained as 


9. (3) Koa Kae es) 
sp eS پچ‎ (1.5) 
9_ (3 SEE в) +5 ` 


A nonlinear ARMA model for the system can be obtained by 
First discretizing the model of Figure I.2 as shown in 


Figure I.3 where F(z) represents the discrete loop filter. 


у (К) 


EN 





Figure I.3. А discrete version of the phase- 
locked loop model. 


-l 


OPENS 





The integration has been approximated as =T 
necessary to use this Euler forward VC SUPR for 
integration rather than one such as the trapezoid rule to 
ayoid a delay free loop). Defining a single linear system 


ag the cascade of F(z) and the discrete integration 


-1 





2 z A(Z) 
КК, ms Brass T-B(z) ul. 5) 
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Equation (4.26) for the phase locked loop becomes 


Жез Ib.OO* o | atot o0 И 

9, (X 0 1 0 0 Жез 

e. = [------------у----------------| |------| (1.7) 
yy G9 0 ШЕ Б БИз SINC) "o 

е(к) -1 | 0 i есю) | 


where it is assumed that 9, (x) and 9, (x) are observed. Using 
the relationships specified by rows 3 and 4, this can be 


written as 


т | Е 
9, (k) b, (x) 0 | 0 0 9, Ck) 
9, (к) 0 ШІ | 0 0 9. (к) 
E = -----------4------------ nn 
Y y °K) 0 0 | 0 Se) Yy tk) 
ЕСЕ) -1 l 0 0 | ES 
a4 Ck) sinl-(.) И o. 20 92 (к) 
| 
+ op En (1.8) 
0 оо у (к) 
0 0 | 


From which it is clear that a finite memory nonlinear ARMA 
model is appropriate, The model can be obtained from the 


first row in equation I.8 by substituting a series expansion 
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for the sin function and truncating it at the degree desired 
resulting in 


M 


С 
(2 


Dot 








9. (k)=b, (k)*9, (Kk) и a, (k)* 2, 


m=0 
ЕС) 
мһеге 2М%1 is the degree of the series approximation to 
sin(.). (Note that lim A, (2) = lim В. (2) = 0 so that the 
right hand side of Г. се О а чае 
ОО Ле output 95 (K).) An infinite memory Volterra series 


model for this system has been considered by Van Trees 


Meee. 55]. 
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Е ВА, J 


Model Simulation Program Listings 


This appendix provides a listing of the fortran model 
Simulation programs used in the experimental study of the 
lattice characteristics, Included are the main programs 
for the batch processing ARMA lattice, the adaptive ARMA 
lattice and the brute force model solution. Each main 
program is followed by a collection of subroutines used only 
by that specific main program. Then a collection of common 
subroutines called from two or more locations in the batch 


lattice, adaptive lattice or brute force method is listed. 
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CALL ПОСРЕУ(К«А,В, ЕК, ВК, АМУР, АЦЕ, АУЕ , АЧВ, ВУЕ, РР, РУВ, ВЕР ,Т, С) 


СЕСЕВ 


ERRCR FCF THIS 
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